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This  report  is  based  on  work  completed  in  partial 
fulfillment  of  the  requirements  for  the  Ph.D.  degree 
in  Biomedical  Fngineering  for  Drexe'l  University 3 
Philadelphia }  Pennsylvania .  Portions  of  the  experi¬ 
mental  results  are  repeated  (for  clarity)  from  the 
Annual  Progress  Report  "Development  of  Multi-Layer 
Models  for  the  Head  under  Impact"  (Contract  #N62269- 
72-C-0171)  November  8}  1972  by  G.D.  Moskowitz  and 
J.L.  Rose. 


SUMMARY 


As  a  consequence  of  Man's  travel  at  increasing  speeds,  the  incidence  and 
severity  of  head  injuries  has  also  increased.  The  development  of  adequate 
protection  against  such  injury  requires  a  thorough  knowledge  of  the  mechanics 
of  trauma  to  an  unprotected  head. 

A  finite  difference  form  of  the  governing  equations  of  motion  for  one  and 
two-dimensional  wave  propagation  is  utilized  to  solve  the  problem  of  non-pene¬ 
trating  impact  to  the  human  head.  Two  aspects  of  the  analyses  provide  signi¬ 
ficant  improvement  over  previous  models:  (1)  multi-layered  structuring  which 
includes  (2)  more  accurate  material  definitions  for  each  layer.  A  series  of 
dynamic  photoelastic  experiments  add  qualitative  confidence  in  the  solution 
techniques . 

The  layered  plate  one-dimensional  analysis  provides  a  method  of  predicting 
the  influence  of  several  material  property  and  size  modifications  in  a  geo¬ 
metrically  simplified  head  impact  model.  The  spherical  model  with  a  layered 
energy  absorbing  skull  yields  highly  attenuated  and  smoothed  tensile  pressure 
peaks  in  the  brain  as  compared  to  the  results  with  a  single  layered  elastic 
skull.  An  elastic  brain  model  (that  includes  an  assumed  high  dynamic  shear 
modulus)  suggests  that  the  combined  shear-normal  stress  levels  would  be  move 
likely  to  cause  failure  than  the  shear  free  stress  condition  in  a  hydrodynamic 
brain  model.  The  generality  of  the  solution  techniques  would  readily  permit 
extension  of  the  analyses  to  investigate  the  significance  of  future  modelling 
considerations . 
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I  -  INTTRODUCTION 

As  a  consequence  of  man's  travel  at  increasing  speeds,  the  incidence 
and  severity  of  head  injuries  has  also  increased.  Walker  (1966)  noted 
that  in  motor  vehicle  accidents  approximately  two-thirds  of  the  cases 
involve  head  injury.  In  addition,  between  the  ages  of  1  and  44  years, 
head  injury  ranks  as  the  prime  killer  in  the  United  States  (Walker  (1966). 
The  development  of  adequate  protection  against  such  injury  requires  a 
thorough  knowledge  of  the  mechanics  of  trauma  to  an  unprotected  head. 

In  till b  present  work,  theoretical  one-and  two-dimensional  head 
models  are  developed  which  could  lead  to  .improved  understanding  of  head 
injury  mechanics.  Two  aspects  of  the  analysis  provide  significant  im- 
p^/onent  over  previous  models:  (1)  multi-layered  structuring  which 
includes  (2)  more  accurate  material  definitions  for  each  layer.  The 
realistic  head  impact  model  solution  predicts  the  stress  wave  propa¬ 
gation  which  occurs  when  the  head  is  struck  by  a  low  velocity  impactor. 

Physiological  and  Pathological  Considerations 

There  are  three  classifications  of  types  of  blows  to  the  head.  Im¬ 
pulsive  blows  involve  accelerations  to  the  head  when  there  is  no  direct 
contact  with  the  impactor.  This  situation  is  typified  by  whiplash 
injuries.  Objects  which  strike  the  head  and  pierce  the  brain  matter 
produce  penetrating  blows.  Bullet  and  high-velocity  fragment  wounds  are 
examples  of  this  type  of  injury.  Finally,  impactors  which  strike  the 
head,  but  do  not  enter  the  brain,  are  classified  as  non-penetrating 
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blows.  Unterhamscheidt.  and  Sellier  (1965)  further  subclassify  non- 
penetrating  blows  as  primary  (cerebral  cortex  damage  and  subdural  hema¬ 
tomas)  and  secondary  (circulatory  problems  due  to  brain  stem  damage) . 

The  model  discussed  in  this  dissertation  is  for  the  case  of  a  primary 
non-penetrating  blow  to  the  head  where  the  skull  does  not  fracture. 

In  the  present  research,  initial  wave  propagation  response  analysis 
is  the  primary  objective;  howevei ,  other  potential  injury  mechanisms  are 
included  to  provide  a  complete  description.  Because  of  the  relatively 
high  compressive  strength  or?  brain  tissue,  overpressure  is  generally 
regarded  as  a  secondary  cause  of  injurious  mechanical  stress.  Potential¬ 
ly  dangerous  forces  resulting  from  underpressure,  rotation,  and  cervical 
stretching  will  be  discussed  individually. 

Benedict  (1969)  discussed  three  means  by  which  underpressure  can 
result  in  tissue  damage:  (1)  outgassing  of  entrained  gas;  (2)  vaporiza¬ 
tion  of  liquid  (flash  boiling);  and  (3)  tissue  failure  in  tension.  Liu, 
Chan,  and  Nelson  (1971)  describe  a  failure  criterion  for  capillary  vessels 
as  th;  cumulative  time  spent  beyond  a  critical  tensile  pressure. 
Unterhamscheidt  and  Sellier  (1966)  showed  that  a  fixed  head  (of  experi¬ 
mental  animals)  allowed  no  separation  between  brain  and  skull,  and 
therefore  produced  no  cavitation  injury.  Goldsmith  (.1966)  expressed  an 
opposing  view,  claiming  that  cavitation  is  the  result  of  rarefaction 
waves  reflected  from  the  skull  boundary  opposite  from  the  site  of  the 
blow. 

Rotational  forces  which  result  in  the  shearing  of  tissues  can  be 
caused  by  whiplash  or  non -penetrating  blows.  The  original  work  in  this 
area  was  presented  by  Holboum  (1943),  who  based  much  of  his  theory  on 
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the  fact  that  the  bulk  modulus  of  brain  tissue  is  several  orders  of 
magnitude  higher  than  the  modulus  of  rigidity.  This  fact  suggested 
that  at  low  shear  levels  brain  tissues  could  separate  and  produce  in¬ 
jury.  Qnmaya,  Hirsch,  and  Martinez  (1966)  subjected  experimental 
animals  to  blows  which  could  cause  head  rotation.  In  one  case,  a 
plaster  collar  restricted  head  rotation  and  permitted  only  translation, 
and  in  the  second  case,  unrestricted  head  rotation  was  allowed.  The 
results  indicated  that  rotation  was  necessary  to  cause  concussion  at  a 
given  level  of  impact.  One  experimental  study  induced  purely  rotational 
forces  on  the  head.  This  resulted  not  only  in  brain  trauma,  but  also  in 
spinal  cord  injury  (Unterhamscheidt  and  Higgins  (1969)). 

Cervical  stretch  has  been  demonstrated  to  produce  injury  patterns 
similar  to  impact  without  any  actual  impacting  blow  (Friede  (1961))  and 
vonGierke  (1966) .  vonGierke  (1966)  used  a  nearly  hardened  soft  plastic 
model  of  the  cervical  cord  region  to  show  the  indentation  caused  by 
stretching  over  the  odontoid  process. 

Injury  produced  by  a  change  in  volume  within  the  brain  case  is 
mainly  related  to  slower  acting  forms  of  insult,  such  as  tumors  or  hema¬ 
tomas  (e.g.,  Evans  (1966)  and  Gurdjian  and  Webster  (1958).  Langfitt 
et  al  (1966)  demonstrated  that  for  monkeys,  an  intracranial  volume 
change  greater  than  4  ml  resulted  in  extremely  high  intercranial 
pressure.  While  this  is  not  of  direct  consequence  in  impact  studies, 
the  cumulative  volume  change  effects  of  repeated  (vibration)  blows  might 
be  important. 

The  neuropatho logical  investigations  of  Lindenberg  and  I-'reytag  (1957) 
show  the  typical  wedge  shaped  contusion  necroses  found  on  the  opposite 
from  the  impact  location.  These  contrc  coup  regions  are  widest  at  the 
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crest  of  a  convolution  with  the  point  of  the  wedge  extending  toward  or 
into  the  white  matter.  Streak-like,  multiple,  and  densely  arranged  con¬ 
tusion  hemorrhages  are  also  described.  Detailed  light  microscopic 
analysis  of  nerve  fiber  (white-matter)  damage  is  given  by  Stritch 
(1956  and  1961).  Evidence  of  tom  or  severed  segments  of  white-matter 
is  shown  with  the  presence  of  retraction  balls  of  axoplasm  at  the  cut 
ends.  A  general  discussion  of  the  feasibility  of  using  piiysiopatho- 
logical  discoveries  to  uncover  the  mechanical  pathways  of  head  injury 
is  presented  by  Oimaya  (1969)  . 

Injury  Tolerance  Criteria 

Experimental  in-vivo  observation  of  neural  tissue  injury  in  humans 
during  head  impact  has  not  been  (and  probably  never  will  be)  accomplished. 
Physical  and  theoretical  modelling  of  head  impact  must  be  cautiously 
correlated  to  the  best  available  in-vitro  and  pathological  human  findings 
and  in-vivo  animal  experiments. 

An  important  consideration  in  modelling  head  impact  is  the  dynamic 
nature  of  a  non-penetrating  blow.  The  duration  of  an  impact  blow  has 
physical  as  well  as  physiological  significance  in  defining  head  injury 
criteria.  Komhauser  (1954)  describes  a  cross-over  point  between  short 
and  long  blows  at  one- third  to  one-fourth  of  the  natural  period  for 
simple  single  degree  of  freedom  systems,  Melbourne  (1943)  and  Rayne 
and  Maslen  (1969)  state  that  long  duration  blows  (compared  to  the 
natural  period  of  the  system)  result  in  head  displacement  being  pro¬ 
portional  to  the  peak  input  acceleration,  while  for  short  duration 
blows  displacement  is  proportional  to  change  in  velocity  . 

The  Wayne-State  Impact  Tolerance  curve  was  developed  with  the 
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criteria  .eventing  a  simple  linear  skull  fracture  that  often 
accompanies  concussion.  Eased  primarily  on  results  determined  by  dropping 
cadavers  on  their  heads,  the  Wayne-State  curves  (e.g.,  Gurdjian,  et  al. 
(1964]  and  Patrick  (1966)  have  the  general  form  of  a  rectangular  hyper¬ 
bola.  The  ordinate  of  this  curve  is  acceleration  and  the  abcissa  is 
duration  of  blow.  Beyond  a  crossover  duration  of  about  10m  sec,  the 
injury  tolerance  can  be  defined  as  a  constant  level  of  peak  acceleration. 

TVo  additional  injury  criteria  are  based  on  the  concepts  of  the 
Wayne-State  curve.  The  Severity  Index  is  computed  by  integrating  (with 
respect  to  time)  acceleration  to  the  2.5  power.  The  single  Tesultant 
value  must  be  less  than  a  given  injury  level  constant  (1000) .  A  recent 
Maximum  Strain  Criteria  by  McElhaney,  Roberts,  and  Stalmaker  (1971) 
describes  the  largest  strain  permissible  in  a  simple  mass -spring -damper 
model  for  head  impact  response.  The  Maximum  Strain  Criteria  is  also 
based  on  the  average  strain  required  to  cause  skull'  fracture. 

One  rotational  injury  criterion  is  not  based  on  preventing  skull 
fracture.  Mahone,  et  al  (1967)  present  rotational  velocity  vs.  rotation¬ 
al  acceleration  injury  threshold  curves  for  experimental  impacts  to 
Rhesus  monkeys.  A  scaling  factor  between  various  primates  was  developed 
by  Chmaya,  Yamell,  Hirsch,  and  Harris  (1967)  such  that  an  inverse  2/3 
power  relationship  exists  between  brain  masses  and  rotational  accelera¬ 
tions.  The  monkey  to  man  brain  mass  similitude  analysis  can  be  used  to 
plot  a  rotational  tolerance  curve  for  man. 

The  Wayne  State  curves  were  developed  by  considering  that  simple 
linear  skull  fracture  frequently  appears  in  concussed  patients.  This 
phenomenological  approach  does  not  account  for  the  mechanisms  of  the 
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primary  concussion  injury.  Likewise,  the  Severity  Index  defines  a 
somewhat  arbitrary  injury  cutoff  value  which  is  not  determined  from 
potential  brain  tissue  damage.  Even  the  rotational  injury  criteria  are 
not  related  to  specific  shear  levels  within  the  cranial  vault.  All  of 
these  criteria  are  useful  in  estimating  a  damage  threshold;  however, 
they  do  not  define  the  severity  or  nature  of  the  brain  injury  mechanism. 

The  Maximum  Strain  Criteria  provides  a  meaningful  correlation  be¬ 
tween  the  existing  protective  indices  and  the  actual  strain  levels  in  a 
lumped  parameter  skull-brain  model.  The  Cumulative  Damage  Index  (liu, 
et  al  (1971))  describes  contre-coup  injury  as  the  time  averaged  stress 
during  which  the  negative  pressure  at  a  given  brain  region  is  above  a 
critical  level.  In  a  like  manner,  the  present  research  effort  is 
directed  towards  defining  physical  phenomena  which  are  the  potential 
source  of  brain  damage  during  impact.  Whenever  possible,  the  brain 
injury  mechanism  hypotheses  are  correlated  to  the  exisiting,  more 
directly  measurable,  impact  tolerance  criteria. 

Time  Scale 

The  term  non -penetrating  head  trauma  suggests  that  at  some  time, 
after  the  blow  is  delivered,  an  injurious  mechanical  force  is  applied 
to  the  brain  for  a  finite  period  of  time.  Two  criteria  are  necessary 
to  define  a  time  scale  for  identifying  types  of  head  injury  mechanisms: 
(1)  time  of  injury  occurence  and  (2)  duration  of  damaging  force.  A 
simplifying  fact  is  that  most  short  duration  injury  forces  occur  soon 
after  the  impact,  while  most  longer  duration  forces  occur  later. 

The  sound  wave  transit  time  in  a  brain  is  about  120  u  sec 
(Unterhamscheidt  and  Sellier  (1966)).  IXiring  the  wave  propagation 
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which  follows  an  impact  blow  many  reflected  and  refracted  waves 
periodically  interact  with  each  other.  This  indicates  that  local 
brain  tissues  are  stressed  with  a  period  several  times  less  than  the 
characteristic  period  of  120  n  sec. 

Only  physical  modelling  experiments  have  shown  the  existence  of 
potential  injury  forces  at  times  less  than  1  m  sec  (1000  u  sec)  after 
impact.  Janssen  and  Bowman  (1970)  show  wave  propagation  in  a  photo- 
elastic.  brain  model.  Because  the  high  shear  stress  levels  enter  and 
exit  a  given  location  from  frame  to  frame  (250w  sec  intervals) ,  the 
damaging  force  duration  must  be  less  than  250  v  sec.  Gross  (1958) 
analyzed  high-speed  photographs  (4000  frames  per  second)  of  a  rubber 
mallet  striking  a  water  filled  glass  flask  to  display  cavitation 
bubbles  500  u  sec  after  impact.  Sellier  and  Unterhamscheidt  (1965) 
implanted  pressure  transducers  in  cadaver  heads  and  recorded  a  peak 
negative  contre-coup  pressure  at  500  w  sec  after  delivery  of  an  impact 
blow. 

Many  of  the  later  occurring  peak  forces  were  observed  by  delivering 
a  blow  to  the  head  of  an  animal  subject  whose  head  was  free  to  rotate. 
Data  from  such  experiments  was  gathered  from  probes  located  external 
to  the  skull.  Friede  (1961)  struck  the  occipital  region  of  a  cat's 
skull  with  a  pendulum  and  recorded  maximal  frontal  bone  velocity  at 
4  m  sec.  Chmaya,  Hirsh,  and  Martinez  (1966)  impacted  monkeys  and  photo¬ 
graphically  determined  the  peak  rotational  acceleration  to  occur  at 
4  msec  and  to  have  a  duration  of  2.5  msec.  In  experiments  with  human 
volunteers  wearing  helmets,  Lombard  (1949)  showed  peak  acceleration 
from  a  frontal  blow  to  last  approximately  2  m  sec. 
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Even  longer  durations  and  later  arrival  times  for  injury  forces  ar^ 
found  in  the  indirect  whiplash  type  of  injury.  Horizontal  sled  experi¬ 
ments  with  humans  (e.g.  Ewing,  Thomas,  Patrick,  Beeler,  Smith  (1970)  and 
Stapp  (1957)  resulted  in  approximately  10  a  sec  duration  peak  head  ac¬ 
celerations.  Head  motion  started  75  m  sec  after  the  onset  of  sled  ac¬ 
celeration. 

The  present  detailed  wave  propagation  analysis  is  not  extended  for 
a  long  time  after  the  impact  blow  is  delivered.  Longer  solutions  to  the 
complex  formulations  of  multi-layered  head  impact  might  require  an  in¬ 
creased  degree  of  approximation.  Emphasis  is  placed  on  predicting  tissue 
damage  phenomena  occurring  earlier  than  500  psec  and  lasting  less  than 
25  usee. 
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II  -  MODELLING  CONSIDERATIONS 
Previous  Theoretical  Models 

Before  beginning  a  thorough  discussion  of  the  present  theoretic¬ 
al  modelling,  a  brief  review  of  previous  skull -brain  impact  analyses 
is  presented.  The  mathematical  approach  of  most  investigators  has 
been  to  consider  geometrically  and  materially  simplified  models  which 
represent  the  overall  closed  head  (no- fracture)  impact  situation. 

This  approach  has  led  to  the  prediction  oi  potential  sites  of  contre- 
coup,  intermediate  coup,  and  coup. 

The  first  theoretical  model  (Anzelius  (1943))  consisted  of  a 
rigid,  fluid-filled  spherical  shell  which  reduced  the  problem  to  the 
solution  of  the  wave  equation  in  spherical  coordinates.  Goldsmith 
(1966)  suggested  that  the  model  be  modified  to  include  an  elastic  shell. 

The  mathematical  contribution  of  Rand  and  DiMaggio  (1967)  in  the 
field  of  acoustics  led  to  a  solution  for  the  vibratory  response  of  a 
fluid-filled,  elastic  membrane  sphere.  Engin  (1969)  added  bending  to 
the  formulation  and  described  the  head  as  a  thin,  elastic,  isotropic, 
and  homogeneous  sphere  filled  with  an  inviscid,  compressible  fluid. 

The  solution  technique  consisted  of  a  Laplace  transformation  method 
with  a  radially  directed  impulsive  loading  function. 

Several  modifications  of  the  Engin  model  have  been  made  within  the 
same  general  configuration.  Benedict  (1969)  considered  the  Engin 
approach  too  cumbersome  for  any  generalized  forcing  function,  except 
the  dirac-delta.  He  has  obtained  a  finite  difference  solution  for  a 
time  and  spatially  varying  input  function.  Recently,  Liu,  et  al  (1971) 
expanded  the  input  to  the  Engin  model  to  include  a  spatially  dependent, 
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finite  time  pulse.  An  approximate  thick  shell  model  with  transverse 
shear  and  rotatory  inertia  was  solved  by  Advani  and  Lee  (1970) .  A 
suggested  improvement  (Engin  ar,d  Liu  (1970)  and  Roberts,  Hodgson,  and 
Thomas  (1966)  is  the  addition  of  a  weakly  plugged  foramen  magnum. 

One-dimensional  models  are  inherently  much  simpler  to  analyze 
and  much  less  expensive  in  computer  time  costs.  The  small  number  of 
parameters  required  to  define  the  problem  also  make  discussion  of  one¬ 
dimensional  models  much  easier. 

Hayashi  (1969)  describes  a  rigid  bu"  massless  container  (skull) 
filled  with  an  elastic  fluid  (brain) .  The  vessel  is  attached  to  a 
linear  spring  which  represents  the  resultant  effects  of  skull  and 
hair  elasticity.  The  slowly  converging  infinite  series  solutions  are 
meaningless  in  computing  numerical  results.  Liu  (1970),  in  reviewing 
and  critiquing  the  Hayashi  model,  provides  an  exact,  closed  form 
solution  for  this  one-dimensional  model.  He  also  suggests  including  a 
damper  in  parallel  with  the  spring.  The  lumped  parameter  head  model  by 
McElhaney,  et.  al  (1970)  provides  linear  brain  displacement  results  for 
various  input  configurations.  The  simplicity  of  this  four  element 
model  (two  masses,  one  spring,  and  one  damper)  allows  solutions  to  be 
extended  to  several  hundred  milli  seconds. 

Anatomical  and  Theoretical  Structuring 

Anatomical  observation  suggests  that  the  head  consists  of  many 
layers  of  different  materials.  A  magnified  view  of  a  head  cross- 
section  would  show  ten  or  more  separate  layers.  A  somewhat  simplified 
head  cross-sectional  sketch  is  shown  in  Figure  1.  Many  layers  are 


included  within  each  of  the  four  regions  indicated.  The  rationale 
for  choosing  the  particular  layers  that  are  included  in  the  theoretic¬ 
al  analysis  is  discussed  in  the  following  paragraphs. 

Scalp 
Skull 

Meninges  and 
Cerebrospinal 
Fluid 

Groin 

Figure  1.  Diagram  of  cross-section  of  a  portion  of  the  head. 

Galford  and  McElhaney  (1970)  claim  the  scalp  acts  as  a  cushion 
in  the  transmission  of  forces  to  the  brain.  Because- the  soft  visco¬ 
elastic  scalp  is  only  partially  restrained,  it  is  free  to  be  pushed 
aside  and  be  readily  displaced  by  a  rapidly  moving  impactor.  The 
scalp  primarily  provides  protection  against  lower  velocity,  glancing 
blows.  While  the  scalp  might  slightly  reduce  the  energy  delivered  to 
the  skill,  the  scalp  is  considered  to  be  of  secondary  importance  com¬ 
pared  to  the  skull.  Therefore,  the  scalp  is  not  included  in  the 
initial  modelling. 

The  strong,  bony  skull  provides  the  most  important  brain  protection 
layer.  The  skull  structure  consists  of  three  layers  the  outer  table, 
the  diploe,  and  the  inner  table.  The  outer  and  inner  tables  are  con¬ 
structed  of  dense,  compact  bone.  The  diploe  region  is  a  weaker,  porous 
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structure  that  contains  many  fluid- filled  interconnected  cavities. 
Inclusion  of  three  skull  layers,  each  with  their  own  properties,  is  a 
meaningful  improvement  over  previous  single -layered  skull  head  impact 
models. 

The  meningeal  region  (dura  mater  is  the  most  significant  sub-  . 
structure)  is  a  physically  tough  and  leathery  layer  which  might 
affect  the  wave  reflections  inside  the  skull.  Also  the  dura  provides 
penetration  protection  to  the  brain  against  the  bony  interior  skull 
ridges.  Qmvaya  (1968)  demonstrated  that  piercing  of  dura  occurred  at 
SOI  of  the  force  level  required  to  penetrate  human  skin.  To  date, 
there  have  been  no  experimental  evaluations  of  the  structural  proper¬ 
ties  necessary  to  define  the  behavior  of  the  meningeal  region. 
Therefore,  it  is  impossible  to  establish  an  approximate  description 
of  this  portion  of  the  head.  Additionally,  the  potentially  complex 
constitutive  relationships  that  might  govern  the  structural  nature 
of  the  meninges  would  add  substantial  difficulties  to  the  mathematical 
formulation.  Presently  no  meningeal  layer  is  included;  however,  the 
need  for  further  evaluation  of  this  region  is  clearly  recognized. 

Head  impact  in  the  absence  of  cerebrospinal  fluid  (CSF)  was  shown 
by  Qmnaya  (1968)  to  result  in  concussion  with  very  low  energy  blows. 
Because  CSF  is  rather  inviscid,  it  can  not  act  as  a  viscous  shock  ab¬ 
sorber.  The  protective  role  of  CSF  seems  to  be  filling  the  portions 
of  the  cranial  vault  not  filled  with  brain  tissue  so  that  the  entire 
internal  cavity  remains  without  void  space.  Engin  and  Liu  (1970) 
discuss  the  similar  densities  of  brain  and  CSF  (1.150  gm/au^  to 
1.007  gm/cm^) ,  In  addition,  the  bulk  modulus,  of  brain  and  CSF  are 
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both  approximately  that  of  water  (2.07x10 


10 


te-). 


Mathematically, 


cm 

CSF  and  brain  can  be  considered  to  be  a  single  continuous  structural 
substance.  Physically,  however,  if  the  theoretical  analysis  predicts 
dangerous  mechanical  forces  at  the  extreme  outer  brain  radius,  then 
the  CSF  and  not  the  brain  will  be  damaged. 


The  theoretical  model  consists  of  the  three  bony  skull  layers  and 
the  brain  region.  A  physical  impactor  (not  an  arbitrary  input  function) 
is  included  to  yield  a  more  realistic  impact  analysis.  Solutions  to 
both  one  and  two-dimensional  formulations  are  presented. 

One -dimensional  analyses  are  typically  simpler  and  less  expensive 
than  two-dimensional  analyses.  The  one -dimensional  head  impact  model 
does  not  include  the  important  closed  container  geometrical  effects  of 
a  spherical  model.  However,  the  wave  propagation  influences  of  a 
layered  structure  with  various  material  properties  can  be  mathematically 
determined.  In  general,  the  one -dimensional  computer  analysis  pro¬ 
vides  an  economical  method  of  investigating  several  parameter  varia¬ 
tions  within  the  basic  layered  head  impact  formulation.  In  addition, 
the  flat  plate  model  can  be  compared  to  continuum  and  lumped  para¬ 
meter  one -dimensional  analyses. 


One-dimensional  particle  motion  requires  eliminating  edge  effects 
by  having  layers  that  are  wide  in  comparison  to  the  thickness.  The 
multi-layered  flyer-plate  class  of  formulations  satisfies  the  general 
problem  description.  Figure  2  shows  the  one -dimensional  impactor, 
skull  layers,  and  brain  region.  All  layer  thicknesses  for  the  head 
are  approximately  equal  to  the  anatomic  values.  Because  exact  one¬ 
dimensional  geometrical  thicknesses  are  not  of  critical  significance. 
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j  a  larger  skull  to  brain  thickness  ratio  was  used  to  provide  a  shorter 

problem  and  more  economical  computer  cost.  The  assumed  sizes  for  im 

} 

i  pactor,  total  skull,  and  brain  layers  are  5.0,  1.5,  and  10.0  cm, 

i 

[ 

>  respectively. 

I 

|  The  cranial  vault  (brain  containing  skull  structure)  appears  to 

i  be  ellipsoidal  in  shape.  Most  coronal  cross-sections  are  basically 

|  circular,  while  the  sagittal  and  transverse  cross-sections  have  the 

t 

|  appearance  of  an  ellipse.  It  is  assumed  that  a  spherical  model  will 

;  provide  a  geometrically  simplified,  but  meaningful  approximation  of 

I  the  head  shape.  The  two-dimensional  layered  head  impact  model  is 

j  shown  in  Figure  3.  The  impactor  is  assumed  to  be  moving  in  an  axial 

direction  along  the  synuietry  axis.  Also  the  head  and  impactor  lie  on 
the  same  axis  of  symmetry.  All  layer  thicknesses  and,  therefore, 
radii  are  appi oxirately  equal  to  the  average  anatomic  values.  The 
frontal  and  occipital  skull  regions  are  used  to  determine  skuli 
thickness  values,  because  these  are  the  most  frequent  sites  of  impact 
injury  and  because  blows  to  the  thinner  temporal  regions  might  pro¬ 
duce  fracture  more  readily.  The  radial  thicknesses  of  the  impactor, 
total  skull,  and  brain  are  5.0,  1.1  and  14.0  an,  respectively. 

Material  Property  Definitions 

A  microscopic  view  of  bone  reveals  long  thin  collagen  fibers 
embedded  in  a  hard  crystalline  matrix  of  hydroxyapatite.  This  de¬ 
scription  suggests  that  bone  materials  behave  as  a  transversely  iso¬ 
tropic  specimen.  Also  a  realistic  analysis  of  brain  tissue  would 
indicate  that  it  is  an  anisotropic  and  inhomogeneous  material. 
Ideally,  one  would  like  to  include  anisiwropv  and  itihomogrneity  in  3 
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head  model;  however,  the  likelihood  of  obtaining  a  solvable  formula¬ 
tion  of  such  a  model  is  small.  Therefore,  the  present  model  assumes 
all  materials  behave  isotropically  and  homogeneously. 

Isotropic  materials  have  three  constants  in  the  constitutive 
(stress -strain)  equations.  Only  two  material  constants  are  inde¬ 
pendent  (Chou  and  Pagano  (1969)).  The  formulations  presented  for 
one-  and  two-dimensional  analysis  require  knowledge  of  the  bulk 
sound  speed  (directly  related  to  bulk  modulus) ,  Poisson’s  ratio, 
and  density. 

The  mechanical  properties  assigned  to  each  of  the  skull  layers 
and  the  brain  are  based,  in  general,  in  the  best  available  experi¬ 
mental  head  materials  evaluation.  In  *968,  Qranaya  surveyed  publica¬ 
tions  on  the  mechanical  properties  of  the  nervous  system.  He  reported 
that  only  three  papers  on  brain  tissue  properties  had  been  presented. 
Also,  very  scant  findings  were  given  for  the  other  "constituents  of 
the  nervous  system.  Under  recent  contracts  supported  by  the  National 
Institute  of  Neurological  Diseases  and  Stroke,  the  West  Virginia  Uni¬ 
versity  Biomechanics  Laboratory,  and  the  Michigan  University  Highway 
Safety  Research  Institute  have  investigated  many  aspects  of  the  material 
behavior  for  scalp,  skull,  dura,  and  brain  tissue.  The  following  dis¬ 
cussion  is  a  brief  sumnary  of  the  pertinent  results  of  these  and  other 
material  property  studies. 

The  inner  and  outer  tables  are  of  similar  construction  and  are 
assigned  the  same  property  values.  Melvin,  Puller,  and  Barodawala 
(1970)  indicate,  in  a  density  curve  for  the  entire  skull,  that  the 
table  density  is  1.94  gm/an  .  While  not  specifically  evaluated  for 
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the  compact  tables,  the  Poisson's  ratio  is  approximated  as  0.3. 

The  West  Virginia  University  investigators  measured  a  sound  wave 
velocity  in  compact  bone  of  3000  meters/sec.  An  elastic  modulus 
of  2.07  x  10  dynes/ aii  at  a  strain  rate  of  100  cm/an/sec  was 

obtained  by  Wood  (1971).  It  is  interesting  to  correlate  the 
material  parameters  of  sound  speed  (CQ) ,  bulk  modulus  (K) ,  density 
(p),  Poisson's  ratio  (  v) ,  and  elastic  modulus  (E).  The  computer 
impact  wave  propagation  analysis  utilizes  the  input  values  Cc= 

3000  meters/sec,  p  =  1.94  gm/an  ,  and  v  *  0.3.  The  following  equa¬ 
tions  are  valid  for  isotropic  materials: 

K  -  pCQ2  (1) 

and  E  *  3K(l-2v) .  (2) 

Substituting  the  input  values  into  equations  (1)  and  (2)  the  elastic 
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modulus  is  computed  as  2.10  X  10  dynes/ cm  ,  which  is  nearly  equal 

to  the  experimentally  determined  value. 

Slight  strain  rate  dependence  has  been  demonstrated  for  hara 
femur  bone  (McElhaney  and  Byars  (1965) .  In  testing  the  human  skull, 
the  workers  at  West  Virginia  University  Biomechanics  Laboratory  suggest 
a  certain  degree  of  rate  sensitivity  in  the  fresher  samples.  Wood 
(1971)  demonstrates  linear  regression  curves  to  show  strain  rate  de¬ 
pendence  of  the  elastic  modulus,  breaking  stress,  and  breaking  strain 
for  compact  cranial  bone.  The  elastic  modulus  increases  by  only  551 
as  the  strain  rate  increases  by  4  orders  of  magnitude  (.01  to  100). 
Therefore,  the  equivalent  of  the  highest  value  for  elastic  modulus 
is  used  as  a  constant,  and  the  material  properties  for  bone  are  not 
considered  to  be  strain  rate  dependent. 
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Random  variations  in  diploe  thickness  do  not  have  a  functional 
relationship  (We st  Virginia  University  Biomechanics  Laboratory  (1970)). 
The  large  standard  deviation  computed  by  McElhaney,  Fogle,  Melvin, 
Haynes,  Roberts,  and  Alan  (1970)  suggests  that  the  use  of  mean  thick¬ 
ness  and  mean  material  property  values  for  the  diploe  might  lead  to 
large  variations  between  impact  sites  and  between  different  individuals. 
The  best  approximation  for  average  properties  was  utilized  in  this 
analysis. 

Melvin  et  al.  (1970)  present  graphic  relations  of  compression 

modulus  versus  diploe  density.  Using  average  density  values,  the 

in  2 

elastic  compressive  modulus  is  1.38  x  10  dynes/ cm  .  The  Poisson's 

ratio  is  reported  as  high  as  0.28  (for  biopsy  specimens  by  West 
Virginia  University  Biomechanics  Laboratory  (1970)  and  as  low  as  0.19 
(by  McElhaney  et  al.  (1970).  An  average  value  of  .25  :s  assumed  for 
the  present  model.  Density  values  for  diploe  have  been  variously  re¬ 
ported  from  1.25  to  2.03  gm/cnC*.  The  lower  value  (Melvin  et  al.  (1970)) 
will  again  be  used  to  maintain  consistency. 

Diploe  sound  speed  measurements  have  resulted  in  paradoxical  con¬ 
clusions.  An  initial  study  by  West  Virginia  University  Biomechanics 
Laboratory  (1970)  and  Martin  and  McElhaney  (1971)  shows  the  diploe 
wave  speed  to  be  faster  than  that  in  the  tables.  The  volumetric 
porosity  of  the  tables  is  about  .38  (large  number  of  tiny  holes)  com¬ 
pared  to  a  value  of  .56  for  the  diploe  (smaller  number  of  large  holes). 
The  West  Virginia  University  Biomechanics  Laboratory  researchers  re¬ 
peated  their  material  evaluation  and  developed  a  porosity  dependent 
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sound  speed  equation.  The  latter  experiments  showed  the  more  physical¬ 
ly  predictable  result  of  a  lower  sound  speed  in  the  diploe  region. 

Using  the  statistical  parameters  suggested  in  the  West  Virginia  model, 
the  sound  speed  is  computed  as  2750  meters/sec. 

The  compressive  failure  strength  for  the  initial  crushing  of  the 
diploe  is  shown  by  Melvin  et  al  (1970)  to  be  only  4.137  x  10  dynes/ cm  , 
This  relatively  low  value  suggests  that  the  diploe  might  crush  or  ex¬ 
hibit  plastic  deformation  during  inpact.  Also  the  crushing  diploe 
struts,  which  form  the  porou-  structure,  can  only  deform  until  the 
cavities  disappear  and  the  diploe  becomes  a  solid  region.  McElhaney 
et  al  (1970)  suggest  a  3  percent  total  strain  over  the  entire  skull 
cross-section  as  a  final  crush-up  strain.  If  the  3  percent  level  is 
modified  to  express  the  fact  that  only  the  diploe  region  is  crushing, 

then  the  final  strain  is  approximately  10  percent.  Therefore,  a  crush- 

9  2 

up  stress  level  is  given  by  1.38  x  10  dynes/ cm  . 

There  is  seme  question  as  to  the  effect  of  strain-rate  on  the  stress- 
strain  relationships  for  the  diploe.  Melvin  et  al  (1970)  indicate  very 
little  scatter  of  the  various  rate -dependent  data  points  in  compressive 
tests  on  diploe.  However,  anatomical  observations  indicate  that  the 
fluid-filled  diploe  cavities  might  result  in  strain-rate  dependence 
in  the  in-vivo  skull.  Because  rate -dependent  effects  are  probably  not 
as  important  as  the  crushing  effects,  a  simpler  rate -independent  diploe 
model  is  used  for  the  initial  analysis. 

Several  mechanical  aspects  of  brain  tissue  have  been  evaluated. 

These  areas  of  investigation  include  parametric  stress -strain  relation¬ 
ships,  bulk  modulus  and  sound  speed  measurements,  and  shear  modulus 
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testing.  An  equation  relating  stress  to  strain  and  strain- rate  was 
determined  from  compressive  tests  performed  by  Estes  and  McElhaney 
(1970) .  Gaiford  and  McElhaney  (1970)  tested  brain  specimens  for 
creep  compliance  and  stress  relaxation.  The  four  lunped  parameters 
of  a  Maxwell -Kelvin  model  are  found  which  best  fit  the  creep  compli¬ 
ance  results.  Driving  point  impedance  characteristics  of  the  head 
were  used  by  Stalmaker,  Fogle,  and  McElhaney  (1970)  to  define  the 
parameters  for  a  mass,  spring,  damper  head  model. 

The  viscous  damper  coefficients  described  in  the  above  investiga¬ 
tions  are  not  appropriate  for  wave  propagation  analysis  require  in¬ 
stantaneous  vs.  equilibrium  moduli!  as  determined  from  the  Hugoniot 
curves.  A  strain-rate  dependent  brain  constitutive  equation  is  in¬ 
corporated  in  tiie  one-dimensional  model,  but  the  necessary  modulii 
are  only  roughly  approximated  in  the  absence  of  experimental  data. 

West  Virginia  University  Biomechanics  Laboratory  (1970)  reported 

3 

an  average  human  brain  density  of  1,05  gm/cm  and  a  bulk  modulus  of 
10  2 

2.10  x  10  dynes/cm  .  Rhesus  monkey  brain  sound  velocity  was  1550 
meter/sec  and  varied  only  slightly  when  measured  across  different 
head  diameters.  Using  the  above  for  K  and  p  in  equation  (1),  the 
sound  velocity  in  human  brain  is  computed  to  be  1450  meter/sec. 

A  statement  by  West  Virginia  University  Bictnechanics  Laboratory  (1970) 
(unsupported  by  tabular  or  graphic  data)  claims  the  dispersion  in 
brain  tissue  to  be  very  small  (S.5t)  up  to  300KKz.  Therefore,  the  bulk 
properties  of  sound  and  modulus  are  essentially  invariant  with  increas¬ 
ing  frequency.  The  1450  meter/sec  sound  speed  for  the  two-dimensional 
model  approximates  the  value  used  in  previous  two-dimensional  head 
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impact  models.  The  bulk  sound  speed  in  the  one -dimensional  model  was 
assumed  to  be  equal  to  the  high  frequency  shock  wave  value  reported 
for  water  in  the  Dynamic  Materials  Property  Library  (Lawrence,  et  al 
(1968)).  This  dynamically  determined  higher  speed  (1893  meter/sec) 
was  utilized  before  the  non-dispersive  brain  statement  was  read.  Be¬ 
cause  of  the  rather  large  expense  in  computer  costs,  the  one-dimension¬ 
al  analysis  was  not  rerun  with  the  lower  sound  speed.  While  the  one¬ 
dimensional  computer  results  arc  not  based  on  a  completely  accurate 
head  model,  the  relative  values  are  still  meaningful  and  valid. 

Fallenstein,  ttjlce,  and  Melvin  (1969)  vibrated  at  (10Hz)  brain 

specimens  between  two  parallel  plates  and  determined  the  shear  storage 

3  2 

modulus  to  lie  between  6-11  x  10  dynes/cm  .  A  shear  storage  modulus 

3  6  2 

increase  from  8.28  x  10  to  1.38  x  10  dynes/an  with  frequency  varia¬ 
tion  from  2  to  400  Hz  was  reported  by  Shuck,  Haynes,  and  Fogle  (1970). 
This  result  was  obtained  by  combining  the  mathematical  solution  for 
torsional  input  to  a  circular  cylinder  and  the  experimental  findings 
of  torque  supplied  to  a  cylindrical  brain  sample.  The  graphic  display 
of  this  data  shows  the  shear  modulus  exponentially  increasing  with 
frequency  (West  Virginia  University  Biomechanics  Laboratory  (197C)). 

A  linear  extrapolation  from  400  Hz  to  the  characteristic  frequency  of 
8300Hz  yields  a  shear  modulus  of  5.50  x  10  dynes/cm  .  The  use  of  an 
exponential  relationship  and  a  higher  frequency  of  brain  tissue  loading 
imply  that  much  higher  values  of  brain  shear  modulus  might  exist  dur¬ 
ing  head  impact. 

The  thread-like  interwoven  nerve  cell  and  nerve  fiber  substructure 
which  comprises  the  brain  might  attain  a  strong  shear  rigidity  under 
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dynamic  loading.  The  shear  modulus  of  weakly  cross -linked  and  dilute 
gel  polymers  increases  by  three  to  five  orders  of  magnitude  with  in¬ 
creasing  frequency  (Ferry  (1961)).  It  is  interesting  to  note  that 
■ost  of  the  polymers  reach  a  peak  shear  plateau  of  1.0  x  10  dynes/cm  . 
If  brain  tissue  performs  in  a  similar  manner,  then  its  dynamic  proper¬ 
ties  must  be  significantly  modified  from  the  low  frequency  properties. 

The  elastic  brain  model  utilized  in  the  two-dimensional  analysis  has  a 

9  ? 

shear  modulus  of  1.0  x  10  dynes/cm  which  results  in  a  Poisson's 
ratio  of  0.475. 

Janssen  and  Bowman  (1970)  used  a  very  soft  urethane  (reported 
tear  strength  is  zero)  as  a  brain  model  for  photoelastic  modelling. 

The  results  showed  stress  waves  propagating  through  the  urethane  be¬ 
fore  the  rotational  effects  could  have  occurred.  If  the  dynamic  shear 
modulus  was  not  higher  tlian  the  low  frequency  value,  then  the  urethane 
could  not  support  shear  and  stress  could  not  be  photoelastically  ob¬ 
served. 

Shear  failure  of  brain  tissue  from  a  rotational  force  is  widely 
accepted  as  an  injury  mechanism  in  the  millisecond  time  regime.  As 
polymeric  materials  change  from  low  to  high  shear  modulus  with  in¬ 
creasing  frequency,  they  also  reach  a  glassy-transition  state  in 
which  they  are  brittle.  A  new  brain  injury  hypothesis  is  proposed 
here:  during  the  initial  impact  response  to  a  direct  blow,  brain 
tissue  can  be  injured  by  a  shear-noimal  fracture  mechanism.  Hus 
injury  concept  also  serves  to  unify  the  influence  of  shear  effects 
throughout  the  entire  time  scale. 
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Experimental  Correlations 

Several  results  of  experimental  head  injury  investigations  have 
been  discussed  previously  in  the  Tune  Scale  and  Injury  Tolerance 
Criteria  sections.  Visual  demonstrations  provide  improved  under¬ 
standing  of  damage  mechanisms.  A  brief  sunmary  of  optical  observa¬ 
tions  of  mechanical  injury  phenomena  is  presented  in  the  following 
paragraphs. 

Optical  analysis  can  be  accomplished  by  photoelastic  and  direct 
observation.  Gross  (1958)  used  high  speed  (4000  frame  per  second) 
photographs  to  show  the  existence  of  regions  of  cavitation  in  water- 
filled  glass  skulls.  Pudenz  and  Sheldon  (1946);  Gosch,  Gooding,  and 
Schneider  (1970) ;  and  Qitnaya  (1966)  have  considered  acute  experiments 
with  animals  which  have  had  lucite  and  lexan  calvaria  implanted. 

After  impact,  the  skull  moves  away  from  the  site  of  the  blow  and  the 
brain  lags  behind  (Orrmaya  (1966)).  The  brain  then' swirls,  while  ro¬ 
tating  about  an  axis  passing  through  the  head's  center  of  gravity. 

The  largest  movements  occur  in  the  fronto-parietal  region. 

The  jjhotoelastic  approach  to  experimental  modelling  of  internal 
stresses  during  head  injury  began  when  Holboum  (1943)  displayed 
pressure  gradients  in  a  gelatin-filled  skull.  Gurdjian  and  Lissner 
(1961)  constructed  a  two-dimensional  head  model.  When  this  model 
was  impacted  from  any  direction,  shear  stress  concentrations  appeared 
at  the  brain  stun.  Flynn  (1966)  used  an  outer  ring  of  high  modulus 
and  inner  disk  of  low  modulus  photo-elastic  material  to  simulate  the 
skull  and  brain,  respectively.  The  recent  work  cf  Janssen  and  Bowman 
(1970)  is  the  most  thorough  attempt  at  photo-elastic  modelling. 
Realistic  sagittal  and  coronal  skull  sections  were  used.  The  location 


24 


and  magnitude  of  the  delivered  blow  was  shown  to  effect  the  resulting 
shear  stress  patterns  in  the  brain.  Also,  the  axis  of  head  rotation 
was  shifted  while  maintaining  a  similar  striking  force.  Finally,  the 
significance  of  dura  attachments  was  demonstrated  by  including,  and 
then  eliminating,  these  attachments  when  molding  the  photoelastic 
model. 

Photoelastic  models  of  a  skull  cross-section  were  constructed  by 
two  research  groups:  Melvin  et  al  (1970)  and  West  Virginia  University 
Biomechanics  Laboratory  (1970) .  These  models  used  a  solid  photo¬ 
elastic  material  to  represent  the  outer  and  inner  tables  and  a  region 
containing  many  holes  to  represent  the  dipioe.  In  both  studies,  a 
static  load  was  applied  to  the  model.  High  stress  concentrations  were 
observed  in  the  dipioe  region,  indicating  this  area  to  be  a  site  of 
potentially  early  failure. 

The  exploding  wire  dynamic  photoelasticity  laboratory  at  Drexel 
University  provides  the  capability  of  observing  wave  propagation  in 
birefringent  models.  An  ultra  high-speed  (lu  sec  between  frame) 
framing  camera  photographs  the  event*  Several  aspects  of  theoretical 
modelling  were  studied  experimentally  to  validate  modelling  assumptions 
and  gain  confidence  in  the  general  applicability  of  the  two-dimensional 
analysis.  Wave  form  and  wave  speed  properties  of  flat  plate  table  and 
dipioe  models  were  investigated.  A  concentric  disk  model  of  a  layered 
skull  and  brain  was  impacted  over  a  polar  cap  to  simulate  the  theoretic¬ 
al  two-dimensional  model. 
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III  -  THEORETICAL  M3DELLING 

The  significant  head  layers  and  their  properties  for  impact  modelling 
have  been  defined  in  the  previous  sections.  Theoretical  one  and  two- 
dimensional  models  require  the  equations  of  motion  for  the  given  geometries 
be  solved  for  the  proper  boundary  and  initial  conditions.  Ihe  geometric 
model  shape  is  divided  into  a  network  of  small  grids.  Finite  difference 
analogs  of  the  governing  partial  differential  equations  are  solved  in 
each  of  the  grid  tones.  The  result  of  satisfying  all  the  equations  in 
all  the  zones  yields  an  approximate  numerical  solution  for  the  head  im¬ 
pact  model. 

Finite  Difference  Formulation  and  Solution 

The  one  and  two-dimensional  solutions  utilize  Sandia  Laboratory  com¬ 
puter  codes  "WQNDY  Ilia"  and  "TOODY  XV”.  The  following  mathematical  de¬ 
scriptions  of  the  formulation  and  computational  procedures  are  based  on 
reports  by  Lawrence  (1970)  and  Benzley,  Berthoff ,  and ‘Clark  (1969).  Be¬ 
cause  the  one -dimensional  analysis  a  degenerate  form  of  the  two- 
dimensional  method,  the  details  of  the  two-dimensional  formulation  are 
described.  Whenever  large  deviations  between  the  two  models  exist,  then 
both  methods  will  be  explained. 

’TOODY  IV"  is  a  Fortran  computer  code  that  solves  the  finite  differ¬ 
ence  analogs  to  the  Langrangian  equations  of  motion  in  two  dimensions. 
Limitation  of  the  model  to  motion  in  two  directions  (essentially  twc  in¬ 
dependent  spatial  variables)  does  not  preclude  the  presence  of  induced 
normal  stress  and  strain  in  a  third  direction.  The  semi-circular,  two- 
dimensional  model  shown  in  Figure  3  is  representative  of  either  a 
cylindrical  or  spherical  cross-section.  The  sphere  of  revolution  is 
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obtained  by  rotating  the  circular  cross-section  about  its  axis  of  symmetry. 
Both  geometrical  configurations  are  shown  in  Figuie  4.  A  computer  input 
parameter  a  determines  the  desired  shape  (6=1  implies  a  cylindrical 
geometry  and  6=2  implies  a  spherical  shape) . 

Shell  or  plate  theories  are  not  incorporated  into  the  formulation. 

The  exact  equations  of  motion  are  solved  by  satisfying  the  laws  of  conser¬ 
vation  of  linear  momentum,  mass,  and  energy.  Also  required  are  the  consti¬ 
tutive  and  strain  rate  -  velocity  equations.  Conservation  of  linear  momen¬ 
tum  in  the  x  and  z  directions  is  expressed  as 
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where  p  is  density 

a  is  acceleration  in  the  given  direction 
is  normal  stress  in  the  x  direction 
normal  stress  in  the  y  direction 
T  is  normal  stress  in  the  z  direction,  and 
T*7'  is  shear  stress  in  the  x-z  plane 
The  conservation  of  mass  equation  is  given  by 

pdv  “  M 


where  M  is  mass  in  an  element,  and 

dv  is  instantaneous  volume  of  an  element 
Energy  conservation  is  described  as 
.  9  E 
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where  F.  is  internal  energy  per  unit  mass 
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P  is  rate  of  work  done  by  pressure  against  a  volume  change,  and 
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Pj  is  rate  of  work  done  by  deviator  stress  against  distortion 
The  strain-rate,  velocity  equations  are  expressed  as 


d**  •  jlh!L 

ax 

d2* « d*2  - 1  (-i“L  ♦  jd L) 

Z  v  a  z  3  x  y 

d22  -  jlhL 

a  z 
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where  d**,  d^,  dzz  is  normal  strain  rate  in  the  given  direction 
d*2  is  shear  strain  rate  in  the  x-z  plane,  and 
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ux,  uz  is  velocity  in  the  giver,  direction 
Equations  (3)  -  (10)  i.-ust  be  satisfied  tor  all  layers  ot  the  models.  An 
additional  equation  is  .equiud,  but  this  constitutive  equation  is  speci¬ 
fied  independently  for  each  individual  la>er.  All  constitutive  relationships 
are  defined  in  the  following  section  entitled  "Constitutive  Equations". 


"T00DY  IV"  defines  the  approximate  spatial  derivative  of  an  arbitrary 
variable  function  (c*)  of  x  and  z  as 
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where  c  it  integration  path  for  line  integral  (^f),  and 
A  is  area  enclosed  by  path. 
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For  example,  using  the  coordinates  in  Figure  5,  the  value  of  a  derivative 


required  for  equation  (3)  or  (4)  (with  c  equal  to  A-B-C-D-A  and  area  ABCDA 
equal  ti  A)  is  given  by 
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where  Tj  is  the  value  of  T  for  quadrilateral  1,  etc. 


Figure  5.  Grid  zone  description  for  finite  difference  algorithm. 

The  value  of  pA  is  obtained  by  averaging  the  values  of  the  neighboring 
regions  to  yield 

pA  =  2  (  ^  +  °2A2  *  c3A3  +  p4A4^  ^  C14) 

where  pj  and  A ^  arc-  the  average  density  and  area  of  quadrilateral  1,  etc. 
Other  areas  or  variable  values  required  are  also  determined  by  linear  inter¬ 
polation  and  averaging  echniques. 
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Time  derivates  are  set  in  finite  difference  form  by  direct  use  of  a 
simple,  second  order,  centered  difference  analog  as  follows: 


,  .  1  n+1/2  ,n+l  ,n 

3  t  J  ~  ~tn+*  -  tn 


(15) 


where  the  superscripts  denote  the  value  of  &  and  time  at  the  prescribed 
time  increments  (n,  n  +  j,  n  +1). 

"WONDY-IIIa"  uses  a  centered  difference  analog  for  spatial  as  well  as 

time  derivatives.  If  the  subscripts  j,  j  +  ,  and  j  +  1  represent  the  end 

£ 

points  and  midpoint  of  the  jth  spatial  zone,  then 


Linear  interpolation  expressions  for  time  and  spatial  zones  are  also  utilized 
in  determining  the  final  difference  equations.  The  complete  finite  difference 
system  of  equations  for  "WONDY-IIIa"  and  "TOODY-IV"  is  presented  in  Appendix  A. 

All  physical  materials  and  boundary  regions  are  defined  by  "TOODY-IV”  as 
projections  on  a  rectangular  index  space.  The  column  (i)  and  row  (j)  identi¬ 
fication  of  each  zone  is  fixed  for  easy  referencing.  The  same  material 
particles  are  always  contained  in  a  given  zone  (a  requirement  of  the  Lagrangian 
system  of  equations) .  There  is  no  limitation  on  the  motion  and  distortion 
allowed  for  any  zone.  Variables  in  mesh  3  (see  Figure  5)  are  stored  with  the 
same  index  as  that  applying  to  point  0.  At  a  given  time  and  location,  the 
variables  are  all  being  solved  for  the  value  at  point  0. 

In  addition  to  the  equations  of  motion*  the  formulation  must  satisfy  the 
proper  boundary  and  initial  conditions.  The  initial  values  for  all  variables 
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are  determined  by  input  card  parameters.  The  boundary  criteria  are  establish¬ 
ed  by  defining  a  one  zone  layer  that  has  the  desired  boundary  effects.  Figure  6 
shows  the  index  space  required  tc  model  the  impactor,  three  layered  skull,  and 
brain.  For  example,  along  the  symmetry  axis  a  symmetry  line  is  defined  which 
in  the  "TOODY  IV"  code  automatically  creates  the  mathematical  conditions  for 
allowing  motion  parallel  to  the  axis  but  no  motion  in  the  normal  direction. 
Analogous  techniques  are  utilized  at  the  outer  radius  of  the  skull  by  includ¬ 
ing  a  free  surface  description.  The  slide  line  allows  different  nodal  points 
to  exist  on  both  sides  of  the  same  physical  line.  Slippage  along  the  slide 
line  is  not  permitted  in  the  brain  region,  but  the  impactor  can  slide  along 
the  skull  surface.  Also  the  boundary  layers  are  defined  as  free  surface  and 
symmetry  (a  symmetrical  axis  line) . 

A  physical  shape  is  assigned  to  the  material  meshes  by  using  a  system  of 
transformation  relationships.  The  impactor,  skull,  and  outer  brain  regions 
(labelled  regions  1,  2,  3  and  4)  are  described  by  the  following  equations 

R  =  ro  '  rl  U-y  C*7) 

Xj,i  *  (wo  ‘  wl  +  Xo  (18) 

Zj,i  =  Rcos  Cwo  *  W1  +  Z0  (19) 

where  r  ,  r, ,  w1 ,  x  ,  z  are  non-dimensional  input  constants  that  are  sped- 
fied  for  each  layer.  With  a  proper  choice  of  input  values,  two  adjacent 
regions  with  a  different  number  of  rows  can  be  equal  in  length.  Effectively 
the  column  (i)  locations  determine  the  radius  and  row  (j)  locations  sweep 
out  the  angle,  Because  the  impactor  and  head  model  are  assumed  to  have  a 
mutual  axis  of  symmetry,  only  one  hemisphere  (a  semi-circle  in  cross-section) 
is  analyzed.  The  results  in  the  opposite  hemisphere  are  directly  available  as 
the  syrmetric  image  of  the  portion  that  is  analyzed.  Figure  (7)  shows  the 
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cross-sectional  view  of  the  layered  hemisphere  model  for  head  impact . 

The  central  regions  in  Figure  7  cannot  be  obtained  from  equations 
(17) -(19).  Figure  8  shows  the  transformation  that  exists  between* the  in¬ 
dex  and  physical  space  for  the  left  side  of  the  central  region.  While 
point  4  in  the  index  space  is  surrounded  by  the  normal  rectangular  zone, 
the  same  point  4  in  the  physical  space  is  part  of  a  triangular  grid.  The 
standard  computations  used  for  all  zones  would  result  in  an  undefined 
solution  at  that  location,  because  the  zero  length  of  one  side  of  the 
zone  would  result  in  division  by  zero.  This  computational  anomaly  was 
resolved  by  using  a  linear  interpolation  based  on  the  neighboring  nodes 
as  the  value  for  the  motion  of  point  4. 


Figure  8.  Transformation  from  Index  to  Physical  Space  for  the  Central  Zones 
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With  all  zones  defined,  the  computation  proceeds  as  follows.  At  time*0, 
the  initial  values  define  all  variables  in  all  i,  j  zones.  The  time  is  then 
incremented  by  the  appropriate  initial  time  step  (At).  At  t-At,  the  new  values 
for  all  quantities  are  determined  in  the  1st  colunn  (i*l),  1st  row  in  1st 
colunn  (j=l)  zone.  Then  the  row  identifier  is  advanced  until  the  entire  1st 
column  has  been  solved  at  the  new  time.  The  2nd  column  is  tl  n  advanced  to 
the  new  time,  etc.,  until  all  meshes  have  been  completed.  The  time  is  now  in¬ 
cremented  by  a  new  time  step,  and  the  procedure  continues  until  the  problem 
is  completed  at  a  predetermined  final  time. 

An  important  aspect  of  finite  difference  methods  for  studying  wave  propa¬ 
gation  has  not  been  discussed.  Due  to  material  non-linearities,  compressive 
waves  tend  to  steepen  into  shock  discontinuities.  A  finite  difference  scheme 
requires  that  variable  values  change  by  small  amounts  within  a  given  zone; 
therefore,  the  infinite  shock  slope  would  yield  invalid  computations.  Extreme¬ 
ly  small  mesh  sizes  could  provide  finite  changes  in  the  shock  vicinity,  but 
the  cost  of  running  even  a  short  problem  would  be  prohibitive.  Artificial 
viscosity  is  introduced  to  insure  that  lower  gradients  exist  and  that  the 
finite  difference  algorithm  are  valid.  It  should  be  noted  that  artificial 

-A* 

viscosity  only  operates  in  the  presence  of  compression  waves,  not  rarefaction 
waves.  Two  forms  of  artificial  pressure  (q)  are  incorporated.  A  quadratic 
viscosity  is  given  by 

-ir-)2  <20> 

where  is  an  assumed  constant  length.  The  quadratic  viscosity  is  extremely 
rate  dependent  and  is  effective  near  shock  gradients.  A  linear  viscosity  which 
controls  small  spurious  oscillations  is  given  in  the  form 

q  *  b2  c  (-ft”) 


(21) 
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where  b2  is  assumed  constant  length,  and 
C  is  sound  speed. 

Deviator  components  of  artificial  stress  are  available  but  not  used  in  the 
present  model.  The  values  computed  from  (20)  and  (21)  are  sunned  to  yield  a 
resultant  total  artificial  pressure.  In  the  linear  momentum  and  energy 
equations  (3) ,  (4)  and  (6) ,  the  viscous  stress  values  are  added  to  the  true 
stress  to  obtain  the  final  system  of  equations. 

The  constants  b^  and  b2  are  modified  to  form  non-dimensional  constants 
Bj^  and  B2-  The  values  of  B^  *  1.2  and  B2  =  .06  were  used  as  suggested  by 
Benzley,  et  al  (1969).  A  sample  problem  (using  "WQNBY-IIIa")  with  low  values 
for  the  artificial  viscosity  constants  was  computed.  The  general  result  was 
that  the  stress  and  energy  output  was  the  same  for  norma]  and  low  constant 
values.  At  points  near  the  compressive  wave  peaks,  the  stress  and  energy 
was  less  than  4%  lower  with  the  normal  values  of  and  B2.  Because  the  re¬ 
sults  were  nearly  equivalent,  the  normal  values  were  used  to  project  the 
difference  formulation  from  invalid  computations  near  a  potential  shock  wave. 

Constitutive  Equations 

The  material  property  survey  (in  the  Material  Property  Definitions  sec¬ 
tion)  suggests  that  several  combinations  of  stress  strain  relations  might 
exist  within  the  multi-layered  head  impact  model.  Because  the  stress  levels 
in  the  skull  region  for  this  analysis  are  below  the  yield  point  for  compact 
bone,  the  inner  and  outer  tables  are  always  represented  as  exhibiting  a 
linear  (elastic)  stress -strain  behavior.  The  collapsing  diploe  region  is 
described  by  two  models:  (1)  a  material  that. has  plastic  yielding  beyond  a 
given  elastic  stress  limit  and  (2)  a  crushable  hydrodynamic  foam.  Most  of 
the  skull  models  utilize  a  water-like  hydrodynamic  brain  representation. 
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However,  two  additional  constitutive  relations  were  considered  (1)  a  strain 
rate  dependent  stress  version  and  (2)  an  elastic  model  which  includes  shear 
stress. 

Several  relationships  are  required  to  determine  the  pressure  and  stress 
components  for  a  particular  material  behavior.  Lee  (1969)  suggests  that  in 
cases  of  high  impactor  velocity  (with  result 5 ng  pressures  above  SO  KBar)  the 
equation  of  state  must  include  an  energy  depend ->t  thermodynamic  term. 
"TOODY-IV"  uses  the  Mie-Gruneisen  equation  to  express  /assure  as  a  funct4' 
of  density  and  energy  in  the  following  form. 

P'PK  *  r(E-Eh)p  (22) 

where  fis  Gruneisen  ratio 

pj,  is  reference  pressure  along  the  Hugoniot,  and 
El,  is  reference  energy  along  the  Hugoniot. 

The  relatively  low  impactor  speed  and  the  dearth  of  energy  dependent 
high  speed  data  on  skull  and  brain  tissue  indicate  that,  an  energy  independent 
(theory  of  elasticity)  analysis  is  the  only  potentially  successful  approach. 

To  show  the  accuracy  of  using  energy  independent  constitutive  equations,  a 
sample  problem  was  analyzed  in  one-dimension  with  an  impactor  striking  a  head 
which  had  a  skull  with  aluminum  properties.  The  energy  dependent  part  of 
the  equation  of  state  is  given  by 

p  ■  <>rE  (23) 

In  the  Scanple  problem  the  experimentally  determined  aluminum  constants  were 
available  from  Lawrence,  et  al  (1968).  For  aluminum  the  Gruneisen  ratio  is 
2.1  (typical  range  for  metals  to  water  is  3.0  to  1.0).  The  results  of  this 
sample  problem  showed  the  energy  deperxient  pressure  term  to  be  less  than  11 
of  the  pressure  component  from  the  energy  independent  elasticity  calculations. 
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The  energy  independent  form  of  the  equation  of  state  is  given  by 

P  "  iCon  (24) 

where  Kq  is  adiabacic  bluk  modulus  at  zero  pressure 
n  is  1  ”P0/p,  end 

Po  is  density  at  zeTO  pressure. 

Equation  (24)  is  directly  analogous  to  the  theory  of  elasticity  relationship 
(see  Chou  and  Pagano  (1969)) 

P  *  K  c  (2E) 

where  e  is  the  volumetric  strain. 

The  following  deviator  equations  must  be  used  to  compute  the  stress 

components 


T**  «  dTxx  . 

P 

(26) 

T*z  „  V2 

(27) 

T22  «  dTZZ 
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T**  -  Tzz 

(29) 

'ivher.2  the  terms  preceded  by  the  superscript  "d”  are  the  deviator  stress  com¬ 
ponents.  In  equations  (26)  -  (29)  the  artificial  pressure  terms  are  added 
to  the  right -hand-side  in  o?der  to  compute  the  stress  components.  If  the 
iraterial  is  hydrodynamic  (zero  deviator  streee)  then,  for  example,  equation 
(26)  yields  that  the  normal  stress  (t**)  is  equal  to  the  negative  of  the 
pressure.  The  water-like  brain  model  is  hydrodynamic  in  behavior. 

For  an  elastic-plastic  material  the  stress  deviators  are  non- zero.  The 
stress  (or  stress -rate)  deviators  are  computed  by  applying  Hooke’s  Law  to 
the  strain  (or  strain- rate)  deviators.  The  strain-rate  deviators  (<*dxx,c*d^r. 
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a 


dd2z,  ddxz)  are  defined  by 

dd“  =-<*"♦!  i  4f- 

(30) 

d«  *  \  1  4f- 

V1 5 

(31) 

(32) 

(33) 

Prom  plasticity  theory  the  deviator  strain-rate  is  divided  into 

part  (V*)  ant^  a  plastic  part  (^d**)  given  (for  example)  as 
“  P 

an  elastic 

■•a** .  -V* . ddf 
e  p 

(34) 

Stress  and  stress -rate  are  related  to  the  elastic  s crain  and  strain- rate 
(respectively)  via  Hooke's  Law.  The  deviator  parts  of  the  stress -rate  are 
described  as. 


A02  zw*z  Vz  -  2<Af 

a£ -  0 

aV2  ,  -  drzz)  zA*2 

at 

AIZ  2W*Z  dTX2  zA22 

___ — -  ♦  *  e 


(35) 

(36) 

(37) 


where  G  is  shear  modulus,  and 

w*z  is  rigid  body  rotation 

The  finite  difference  form  of  equations  (35)  -  (37)  can  be  solved  explicitly 
for  the  stress  deviators  at  time  n+1. 

If  the  material  exhibits  plastic  behavior,  then  there  is  an  upper  limit 
to  the  deviator  stresses  whi<~h  is  determined  by  the  yield  condition.  A  typical 
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yield  criterion  is  the  well-kncwn  von  Mises  relationship  (Mendelson  (1968)) 

l£(T3a-Tyy)2+(Tyy-TZZ)2+CTZZ“TXX)2+6(rXy2+Tyz2+TZx2))s  1  y2  (38) 

where  y  is  the  flow  stress  in  uniaxial  stress*  Transforming  this  relationship 
into  deviator  stress  terns  and  substitution  of  equation  (30)  for  yields 

2{(dTxx)2  +  (V2)2  +  (dTZZ)2+(dTXX)(dTZZ)>-  |  y2  (39) 

When  plastic  strain  occurs,  it  is  required  to  be  noimal  to  the  yield  surface. 

The  yield  stress  (y)  may  vary  due  to  strain  hardening,  strain-rate,  or  cnnw'ression. 
However,  in  the  present  diploe  modelling,  y  is  assumed  to  be  constant,  and 
the  material  behavior  is  termed  elastic-perfectiy  plastic. 

The  outer  and  inner  tables  develop  deviator  stresses  but  they  do  not 
exhibit  plasticity.  This  perfectly  elastic  condition  is  accomplished  by 
assuming  an  effectively  infinite  yield  flow  stress.  One  diploe  and  one  brain 
model  also  utilize  a  perfectly  elastic  formulation.  The  former  case  produces 
a  skull  model  that  behaves  as  a  single  linear  elastic  shell  which  serves  as  a 
baseline  for  comparison  with  other  diploe  models  studied  here  and  previous 
investigations  by  others. 

The  collapsing  diploe  cavities  can  also  be  modelled  as  a  crushable.  foam. 

The  equation  of  state  for  hydrogynamic  foams  or  distended  solid  materials  is 
described  by  Herrmann  (1968)  and  Lawrence  (1970).  A  parameter,  a,  describ¬ 
ing  the  degree  of  distention  of  the  foam  is  defined  as 

Ps 

«  “  ~ -  (40) 


where  ps  is  solid  material  density,  and 
p  is  distended  material  density. 
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The  a  parameter  must  vary  to  describe  three  separate  material  regions: 
(1)  initial  elastic  distended  region;  (2)  crushing  distended  region;  and 
(3)  final  elastic  solid  region  (after  final  crushing) .  At  the  boundaries 
between  these  regions  certain  continuity  and  smoothness  conditions  are  im¬ 
posed  on  the  a  parameter.  During  the  crushing,  the  a  parameter  is  assumed 
to  vary  quadratically  with  pressure.  As  in  the  case  of  the  elastic-plastic 
constitutive  equations,  the  pressure  as  well  as  the  a  parameter  are  chosen 
to  be  energy  independent  functions. 

The  rate  dependent  brain  layer  model  is  described  by  a  visco  plastic 
constitutive  equation  for  one-dimensional  wave  propagation  as  described  by 
Hermann,  Lawrence,  and  Mason  (1970).  The  constitutive  relation  is 

-  2Gg  (41) 


where  a  is  total  stress 

oi  is  instantaneous  stress,  and 
g  is  relaxation  rate 


TV«j  of  the  terms  in  equation  (41)  are  further  defined  by 


dqi  _ 


Tt” 


=  Fd 


where  F  is  local  instantaneous  modulus 


d  is  volumetric  strain  rate. 
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Where  o1  is  total  stress  deviator 

a-  £  is  equilibrium  stress  deviator,  and 
tr  is  relaxation  time 

All  required  moduli  are  computed  within  the  progTaa  code  from  the  input  con¬ 
stants  of  density,  bulk  sound  speed,  and  Poisson's  ratio.  If  the  relaxation 
rate  is  small,  then  there  is  not  enough  time  for  plastic  stresses  to  develop. 

On  the  other  hand,  for  large  relaxation  rates  plasticity  occurs  so  rapidly 
that  there  is  not  sufficient  time  for  viscous  over-stress  to  develop.  In  tie 
absence  of  brain  tissue  experimental  data,  the  relaxation  time  is  approximated 
as  one-tenth  of  the  characteristic  wave  transit  time  in  the  brain.  The  one- 
tenth  assumption  is  a  physically  meaningful  value  which  will  show  the  potential 
influence  of  viscous  effects  in  the  brain. 

Stability  and  Accuracy 

The  "WOMDY  III -a"  and  'TOODY  IV”  computer  codes  provide  error  checks  as 
well  as  a  stability  criterion.  Herrmann,  Holzhauser,  and  Thompson  (196/)  de¬ 
veloped  an  approximate  analytic  method  for  determining  a  stable  time  step 
increment  for  equations  of  motion  in  uniaxial  strain.  Valid  time  increments 
are  restricted  by  the  following  stability  equations  in  "WCNDY  Ill-a” 


where  Bj  and  are  artificial  viscosity  constants.  A  similar  stability 
analysis  has  been  extended  to  the  two-dimensional  "TOODY  IV”  rode.  Benzley, 
et  al  (1969)  state  that  this  criterion  has  yielded  efficient  and  stable 
time  steps  for  all  problems  heretofore  solved. 


Iapzctors  cause  large  initial  velocity  gradients;  therefore,  a  smaller 
At  than  required  by  equation  (43)  is  used  at  the  start  of  the  calculation. 

As  problem  time  progresses,  the  value  of  At  is  allowed  to  return  to  the  in¬ 
crement  computed  by  the  stability  equations. 

Using  velocity  interferometer  techniques,  Lundergan  and  Drumheller  (1971) 
achieved  good  correlation  between  experimental  layered  ulate  impact  and  ana¬ 
logous  computations  using  WONDY  Ill.  Experimental  verification  of  a  TOODY  III 
mathematical  analysis  is  presented  by  Karnes  and  Bertholf  (1971).  Other 
experimental  comparisons  to  these  computer  cedes  also  indicate  the  general 
acceptability  of  the  mathematical  analyses. 

The  total  internal  and  kinetic  energy  is  computed  for  all  meshes  at  the 
end  of  each  cycle.  A  check  is  made  to  determine  if  total  energy  has  been 
conserved  while  solving  the  equations  of  motion.  An  instability  or  accuracy 
error  would  result  in  an  energy  blow  up  or  a  stable  but  excessive  error,  re¬ 
spectively.  "WONDY  Ill-a"  also  checks  momentum  conservation  by  a  similar 
procedure.  Unfortunately  the  kinetic  energy  computation  (only  used  in  the 
error  check)  for  the  two-dimensional  analysis  is  the  least  accurate  approxi¬ 
mation  of  the  entire  solution.  For  the  complex  geometry  used,  the  "TOODY  III" 
error  check  is  a  good  measure  of  stability  but  not  a  good  accuracy  checFc  Com- 
pc  ;d  errors  of  up  to  12%  were  considered  acceptable. 

Another  widely  used  method  to  verify  that  stability  and  accuracy  have 
been  achieved  is  to  run  the  same  finite  difference  analog  with  different  grid 
rone  sires.  If  the  two  solutions  are  approximately  equivalent  (within  a  tol¬ 
erable  percentage  difference),  then  the  grid  size  is  assumed  to  be  acceptable. 
The  head  impact  models  were  run  with  zone  dimensions  approximately  one-third 
to  one-half  the  original  zone  size.  The  output,  from  the  smaller  zone  sizes 
produced  results  which  differed  only  slightly  (less  than  5%)  from  the  largei 
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zone  sizes.  Therefore,  the  present  choice  of  zone  sizes  provides  a  valid 
grid  network  for  the  finite  difference  formulation. 


Non-Dimensionalization  and  Input  Value  Summary 

"WQNDY  Ilia"  and  "TOODY  IV"  contain  no  internal  dimensional  constants; 
therefore,  any  self  consistent  set  of  units  can  be  used.  A  properly  non- 
dimensionalized  system  can  be  developed  with  zespect  to  a  characteristic 
velocity  (ur),  density  (p  ) ,  and  length  (xr) .  The  resulting  non-dimensional 
independent  and  dependent  variables  are  given  by 
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The  present  non-dimensional  scheme  is  developed  by  considering  the  size 
and  properties  of  the  brain  (in  the  one-dimensional  analysis)  as  the  character¬ 
istic  values.  These  constants  (based  on  the  one- dimensional  brain  thickness, 

3 

density  and  bulk  sound  speed)  are  given  as  -  10cm,  py  =  l.Ogu/cm  ,  and 
5 

u  =  1.893x10  cm/sec. 
r 
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The  non-dimensional  layer  sizes  are  summarized  in  Table  1.  Because  of 
zone  size  considerations  in  the  two-dimensional  analysis,  the  sirigle  layered 
elastic  skull  has  a  dimensionless  thickness  of  0.10,  while  the  total  thick¬ 
ness  for  the  layered  models  is  O.li.  The  impactor  velocity  is  approximated 
by  the  free  fall  velocity  of  an  object  falling  from  a  high  shelf.  This 
velocity  is  9.465  x  10  cm/sec  or  0.00S  non-dimensional  units.  All  other 
layer  velocities  are  initially  zero. 


One  previous  spherical  head  impact  analysis  has  considered  a  time  and 

spatially  dependent  pressure  input  function.  Benedict  C1969)  defines  the 

time  related  portion  of  the  pressure  (p(t))  as 

2 

P(t)  «  2300  (e'4,73  1  /10'6}(sin(*/10'3)  t)  (45) 

The  force  applied  to  the  head  can  be  directly  determined  as  the  product  of  the 
pressure  times  the  surface  area  of  a  30°  polar  cap  on  the  skull.  The  integral 
of  this  force-time  curve  (for  an  assumed  one  millisecond  duration)  is  equal  to 
the  change  of  momentum  of  the  impactor.  Using  the  impactor  velocity  of  9.465 
x  10  cm/sec,  the  equivalent  mass  required  produce  the  given  force  is 
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computed  as  2520  gms.  Assuming  that  an  iron  impactor  (with  dimensions  de~  i 

l 

scribed  before)  is  utilized  in  the  mathematical  two-dimensional  model,  the  i 

impactor  mass  is  calculated  as  2290  gms.  Therefore,  the  input  energy  is 

i 

approximately  equivalent  for  the  pressure  function  and  impactor  models.  ' 

While  the  force-time  areas  might  be  similar  for  the  two  input  models,  the  ’ 

actual  force  duration  and  rise  time  is  shorter  for  the  impactor  problem. 

i 

The  one-millisecond  duration  pressure  pulse  is  only  attainable  by  consider-  ! 

( 

ing  the  physical  effects  of  a  soft  impactor  (e.g.  a  padded  dash).  The  im- 
pactor  model  is  a  realistic  descriptor  for  a  blow  delivered  to  the  head  model  ; 

presented  in  this  work.  Possible  modifications  to  the  impactor  or  head  model 
will  be  discussed  in  Chapter  VI  -  Recommendations. 

The  present  impactor  model  is  more  closely  related  to  the  dirac-delta 
impulsive  function  utilized  by  Engin  (1969),  and  Liu,  et  al  (1971).  The  iron 
impactor  material  properties  are  those  given  by  Lawrence,  et  al  (1968). 

For  convenient  reference  a  summary  of  the  important  non-dimensional 
material  properi ies  for  the  various  head  layers  used  in  the  two-dimensional 
analysis  are  presented  in  Table  2.  The  elastic  yield  or  initial  crushing 
value  was  decreased  to  one  tenth  of  the  experimental  value.  With  the  orig¬ 
inal  value,  relatively  small  amounts  of  crushing  occurred.  A  larger  impact 
energy  would  induce  more  crushing.  Because  the  experimental  value  might 
be  lower  in  a  dynamic  test,  it  was  more  convenient  to  lower  the  initial 
crushing  value  and  use  the  same  impact  parameters.  The  overall  effect  ir 
to  emphasize  the  energy  absorbing  protection  provided  by  the  skull.  In 
the  one-dimensional  model  the  crush-up  strength  was  inadvertently  set  to  a 
value  of  1.0  which  gives  a  more  gradual,  less  energy  absorbing  pressure- 
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distension  curve  chan  for  the  two-dimensional  case.  Also  the  flat  plate 
models  utilize  different  brain  properties  which  were  described  previously. 


TABLE  2  -  Summary  of  Non-Dimensional  Material  Properties  for 
Head  Impact  Analysis 


NON-DIMENSIONAL 

MATERIAL 

PROPERTY 

1  “  layer  Description  1 

tables 

and 

elastic 

diploe 

elastic- 

plastic 

diploe 

crushable 

diploe 

hydro 

dynamic 

brain 

elastic  | 

brain  1 

density 

1.94 

1,25 

1.25 

1.05 

1.0S 

bulk 

sound  speed 

1.60 

1.467 

1.467 

.766 

.766 

Poisson's  ratio 

0.30 

0.25 

-- 

.475 

elastic  yield 

00 

.00115 

.00115 

OD 

crushup 

strength  j 

-- 

... 

0385 

-- 

-- 

f 

I 

5 

\ 

\ 


fti 
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IV  -  RESULTS  AND  OBSERVATIONS 

Theoretical  Results 

Commuted  results  of  the  one  and  two-dimensional  analyses  of  head  impact 
are  presented  in  condensed  tabular  and  graphic  form.  The  previously  de¬ 
scribed  non-dimensional  units  are  used  for  all  output  variables.  It  should 
be  iiotod  that  the  results  are  based  on  the  impact  response  with  the  assumed 
impactor  size  and  velocity.  Because  of  the  high  computer  costs,  no  attempt 
has  been  made  to  determine  if  the  general  response  characteristics  of  the 
skull -brain  system  are  strongly  dependent  on  the  impactor  parameters. 

Several  computer  output  dependent  variables  such  as  particle  velocity, 
normal  stress,  shear  stress,  wave  speed,  and  energy  density  at  each  spatial 
zone  are  available  for  each  output  time.  These  values  can  be  plotted  as 
variable  versus  position  profiles  at  fixed  times.  Also  a  particular  loca¬ 
tion  can  be  used  to  obtain  a  variable  versus  time  curve.  The  peak  values 
and  duration  above  a  given  cut-off  level  can  be  determined  from  these  out¬ 
put  displays. 

When  the  one-dimensional  analyses  were  performed,  the  computer  soft¬ 
ware  for  automatic  plotting  was  not  available.  Therefore,  the  one -dimension¬ 
al  output  (dependent)  variables  were  printed  in  both  a  complete  and  reduced 
form.  The  reducer!  output  consisted  of  observing  only  seven  locations  within 
the  brain  region,  but  these  new  variable  values  were  sampled  every  0.1  non- 
dimensional  time  units.  The  complete  output  scheme  printed  the  output  vari¬ 
ables  for  every  spatial  zone  in  the  brain  at  evt  y  2.0  units  of  time.  With 
this  latter  output  an  accurate  stress  profile  can  be  drawn  from  the  complete 
spatial  information.  A  sample  stress  profile  in  the  brain  is  shown  in 
Figure  S.  This  curve  shows  the  positive  tensile  stress  levels  at  every  point 


Norm*!  Stress  (non-dimensional  units) 
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in  the  brain  layer  at  time  equals  4.0  non-dimensional  units.  Tables  of 
peak  tensile  stress  values  for  each  one-dimensional  model  are  presented  in 
Chapter  V  -  Discussions. 

A  California  Computer  Flotting  System  was  available  when  the  two- 
dimensional  analyses  were  performed .  Detailed  stress  profiles  and  variables 
versus  time  curves  are  given  in  Chapter  V  -  Discussions. 

A  single  quantity  which  defines  at  each  instant  the  overall  brain  re¬ 
sponse  to  an  impact  would  provide  a  simple  method  of  comparing  the  results 
of  several  models  and  serve  as  a  link  between  continuum  and  lumped-para¬ 
meter  models.  The  internal  energy  in  the  entire  brain  region  yields  a 
measure  of  the  stress-strain  energy  that  has  been  transmitted  to  the  brain. 
The  total  enei’gy  which  equals  internal  energy  plus  kinetic  energy  was 
utilized  for  a  portion  of  the  one-dimensional  analysis.  Specific  curves 
are  presented  in  Chapter  V  -  Discussions. 

Experimental  Observations 

The  analytic  model  provides  the  detailed  mechanical  results  of  a 
theoretical  head  impact,  but  there  is  a  lack  of  experimental  in-vivo 
evidence  to  correlate  these  mathematical  computations  to  physiological 
events.  The  correlation  of  brain  trauma  to  mechanics  is  not  attempted  in 
this  present  work;  however,  experimental  models  are  developed  which  provide 
qualitative  confidence  in  the  analytic  solution.  Dynamic  evaluation  of  two 
aspects  of  head  impact  phenomena  have  been  considered:  (1)  wave  propagation 
through  a  porous,  three-layered  skull,  and  (2)  wave  propagation  in  a  simpli¬ 
fied  concentric  ring  model  for  head  inpact. 

All  experiments  were  performed  at  the  Wave  Propagation  Research  Center 


of  Drexel  University  under  the  direct  supervision  of  Dr.  Joseph  I..  Rose. 
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The  wave  propagation  experiments  utilized  the  techniques  of  dynamic  photo¬ 
elasticity.  A  brief  sisimary  of  the  test  methods  is  discussed  in  the 
following  paragraphs.  A  complete  description  of  the  test  facilities  can 
be  found  in  papers  by  Rose  and  Chou  (1970  and  1972).  The  photographic 
results  are  displayed  in  Giapter  V  -  Discussions. 

Pressure  from  an  exploding  wire  can  be  used  to  impulsively  load  a 
thin  plate  photoelastic  model  of  the  layered  skull,  or  of  the  entire  head. 

The  propagation  of  stress  waves  is  photographed  by  a  Beckman  and  Whitley 
Model  201  framing  camera.  The  camera  aperture  was  reduced  in  size  to 
yield  an  exposure  time  of  approximately  .3u  sec.  This  short  exposure 
minimized  the  blurring  due  to  wave  motion  during  the  exposure  period.  The 
time  duration  between  frames  for  the  present  testing  is  lu  sec. 

Light  is  made  available  from  an  E.G.  and  C.  FX.-1C-6  xenon  flash  tube. 

The  xenon  flash  tube  is  mounted  along  the  focal  line  of  a  cylindrical 
parabolic,  chromium-plated  reflector.  Light  from  the  flash  tube  and  re¬ 
flector  then  enters,  in  order,  a  polarizer,  1/4 -wave  plate,  test  specimen, 

1/4 -wave  plate,  analyze*,  filter  and  camera. 

The  detailed  model  of  the  pcrous  diploe  layer  of  the  skull  is  based  on 
a  photomicrograph  by  Melvin,  et  al  (1970) .  Even  the  hard  table  regions  have 
microscopic  fluid-filled  cavities.  The  effect  of  these  tiny  holes  is 
modelled  by  several  rows  of  small  (.2  an)  drilled  circular  holes.  Figure  10 
shows  a  sketch  of  the  two  skull  models  (the  second  model  was  constructed 
with  more  diploe  details  to  provide  a  lesser  free  path  for  wave  propagation) 
and  the  reference  material  plates.  Hie  reference  section  provides  measure¬ 
ments  of  undisturbed  wave  speeds  for  the  solid  and  void  filler.  The  disk 
and  concentric  ring  head  model  is  shown  in  Figure  11.  Flynn  (1966)  constructed 


b.  Solithane  113  Filler 


Figure  10.  Photoelastic  Skull  Model 
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a  similar  model  which  broke  after  a  few  impacts.  Therefore,  the  skull 
portion  of  the  present  head  model  was  machined  thicker  than  necessary 
for  scaling  to  the  anatomic  size.  Also  a  separate  support  bracket  was 
supplied  to  prevent  the  specimen  from  sagging  under  its  own  weight. 

The  photoelastic  model  mist  have  material  properties  which  are 
similar  to  the  corresponding  regions  of  the  head.  In  addition,  the 
birefringent  nature  of  the  plastic  specimens  must  produce  readily  observ¬ 
able  stress -optical  fringe  patterns.  Wave  speed  (plate  velocity  in  these 
experiments)  and  acoustic  impedance  (wave  speed  x  density)  are  the  critical 
properties  in  wave  propagation  analyses.  In  the  layered  skull  model  the 
acoustic  impedance  at  each  solid-filler  interface  determines  the  propor¬ 
tions  of  reflected  vs.  refracted  wave  propagation.  In  the  head  model  the 
relative  wave  speeds  in  the  outer  ring  and  central  disk,  also  strongly  in¬ 
fluence  the  resultant  wave  motion.  Several  materials  were  utilized  be¬ 
cause  of  their  material  characteristics  and  their  availability  in  the 
required  solid  sheet  or  castable  form.  Table  3  shows  a  comparison  of  the 
important  wave  propagation  properties  for  the  anatomic  and  plastic  model 
materials . 

The  experimental  values  for  the  wave  speeds  in  the  photoelastic 
materials  are  determined  by  measuring  tlte  plate  velocity  for  each  material 
in  the  reference  sheet.  Because  no  two  explosions  are  exactly  the  same, 
all  wave  speeds  must  be  determined  from  the  12  photographs  of  a  single 
test.  By  utilizing  isotropic  elastic  constants,  the  bulk  speed  (CQ)  may  be 
computed.  There  is  a  discrepancy  between  the  static  property  prediction  of 
the  bulk  speed  and  the  dynamic  experimental  values.  The  lack  of  knowledge 
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of  the  wave  speeds  before  testing  makes  the  choice  of  materials  very  diffi¬ 
cult.  In  the  present  testing,  acoustic  impedance  mismatches  were  higher  and 
lower  than  those  for  the  actual  anatomic  regions. 
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V  -  DISCUSSION 

One-Dimensional  Analysis 

The  flat  plate  one-dimensional  modelling  of  the  skull  is  not  directly 
comparable  with  closed  container  spherical  models.  In  addition  to  the  un¬ 
connected  skull  regions,  the  plate  impactor  strikes  the  entire  skull  on 
one  side  of  the  model.  The  energy  density  for  such  an  impact  is  not  the 
same  as  the  case  of  a  small  impactor  hitting  only  a  portion  of  the  skull. 

This  concept  does  provide  a  model  of  wave  propagation  emanating  from  one 
side  of  the  head  and  travelling  back  and  forth  through  the  brain.  The 
primary  objective  of  the  plate  model  is  to  economically  (in  terms  of  com¬ 
puter  expenses)  determine  the  effect  of  several  parameter  variations  with¬ 
in  the  basic  energy  absorbing  multi-layered  skull-brain  system.  Therefore, 
the  results  will  emphasize  the  differences  in  wave  propagation  response  for 
several  types  of  flat  plate  head  models. 

The  nature  of  the  mathematical  coupling  between  the  brain  and  innei 
tables  was  found  to  be  very  important.  Previous  mathematical  models  have 
considered  a  strong  bond  at  the  hrain-inner  table  boundary,  which  caused 
the  particles  on  both  sides  of  the  boundary  to  move  together  without  sepaiat- 
ing.  This  non-separating  connection  permits  the  dissipative  and  dispersive 
effects  of  the  diplce  to  remain  influential  during  the  entire  period  of 
analysis.  A  separable  situation  can  be  modelled  by  assuming  that,  the  mathe¬ 
matical  coupling  between  these  two  layers  can  fail  above  a  given  tensile 

limit.  To  show  this  effect  quantitatively,  an  arbitrary  separation  tensile 
8  2 

stiess  of  1.0  x  10  dynes/cm  was  assumed  for  both  brain-inner  table  boundar¬ 
ies  (near  the  site  of  the  blow  and  on  the  opposite  side).  Separation  occurred 
first  at  the  opposite  side  from  the  blow  (1.68  non-dimensional  time  units) 
and  then  near  the  site  of  the  blow  (2.34  units).  The  separation  (gap)  between 
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the  brain  and  inner  table*  produces  a  high  acoustic  impedance  mismatch  at 
the  brain  boundary;  therefore,  the  stress  waves  are  totally  reflected  back 
into  the  brain  and  not  transmitted  to  the  table  region. 

Figure  12  shows  the  total  energy  (internal  energy  plus  kinetic  energy) 
in  the  brain  as  a  function  of  time  for  the  case  when  the  inner  table  separ- 
ated  from  the  brain  at  a  tensile  stress  of  1.0  x  10  dynes/an  .  Displayed 
in  Figure  13  is  the  total  energy  curve  for  the  case  when  mathematical  separa¬ 
tion  is  not  allowed  (effectively  an  infinite  tensile  stress  for  separation) . 
In  both  Figure  12  and  Figure  13,  the  following  material  properties  for  the 
skull  layers  are  indicated:  elastic  tables  ancl  elastic  diploe;  elastic 
tables  and  elastic -plastic  diploe;  and  elastic  tables  and  crusha'ble  hydro- 
dynamic  foam  diploe.  The  lower  total  energy  levels  (also  lower  internal 
energy  levels),  in  the  case  of  no  separation,  is  due  to  the  continued 
attachment  die  skull  to  the  brain. 

The  complex  C.S.F.  and  meningeal  regions  physically  comprise  what  has 
been  described  as  simply  a  boundary  condition.  Destruction  of  the  vascular 
tissue  in  this  region  could  cause  serious  intracranial  bleeding.  No  experi¬ 
mental  investigations  have  been  performed  have  to  determine  the  nature  or 
level  of  such  injury  mechanisms .  Because  the  one -dimensional  plate  model 
does  not  have  the  geometrical  dispersion  effects  of  a  spherical  model 
(small  impactor  striking  a  large  head) ,  the  stress  levels  are  necessarily 
higher  than  in  a  real  head  impact.  Therefore,  in  the  spherical  models  the 
'tensile  stress  level  at  the  brain- inner  table  boundary  will  be  less  than  in 
the  plate  models,  Until  a  more  definitive  separation  stress  can  be  ascer¬ 
tained,  the  firm  bonding  condition  will  be  utilized  in  all  model  analyses. 
The  potential  importance  of  properly  defining  this  boundary  or  correctly 


incorporating  these  anatomic  layers  requires  future  experimental  and  analytic 
investigation. 

The  internal  energy  within  the  brain  quantifies  the  stress-strain 
tissue  deformation  in  the  brain.  In  comparing  the  effects  of  various  materi¬ 
al  descriptions  for  the  layered  models ,  the  internal  energy  versus  time 
curves  will  be  displayed.  The  curves  in  Figures  14  to  18  are  assumed  to 
have  no  separation.  Because  it  was  desired  to  accurately  plot  data  points 
over  a  large  range  of  values,  a  semi- logarithmic  display  form  was  utilized. 
The  non-dimensional  units,  as  described  previously,  are  used  on  all  curves. 

Figure  14  shows  the  internal  energy  curve  for  the  average  material 
properties  as  listed  in  Table  I.  The  energy  contained  within  the  brain  is 
in  the  following  decreasing  order:  elastic  diploe,  elastic-plastic  diploe, 
and  crushable  foam  diploe.  The  percentage  difference  between  the  elastic 
and  elastic-plastic  models  varies  from  2%  to  50%,  while  there  is  a  10%  to 
65%  difference  between  the  elastic  and  crushable  models.  The  crushable 
foam  model  would  absorb  more  energy  if  the  final  crush-up  pressure  were 
not  set  higher  than  the  actual  value.  Tue  periodic  peaking  of  the  internal 
energy  represents  the  movement  of  the  main  wave  front  as  it  carries  a  large 
portion  of  the  energy  from  the  brain  compartment  to  the  skull  and  back  again. 

The  thickness  of  the  diploe  in  proportion  to  the  table  thickness  differs 
between  locations  on  the  skull  and  between  various  individuals.  A  special 
situation  with  the  diploe  thickness  increased  by  2-1/2  times  the  original 
value  (0.12S  non-dimensional  units  from  0.0S  units)  was  analytically  examined 
to  determine  the  effects  of  a  large  change  in  the  diploe  thickness  parameter. 
Figure  (15)  shows  the  results  for  the  elastic-plastic  and  crushable  foam  en¬ 
larged  diploe  models.  The  elastic  diploe  model  results  are  included  for 
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reference  to  the  original  problem  values.  The  percentage  difference  between 
the  elastic  and  elastic-plastic  models  is  now  as  high  as  651,  while  the  per¬ 
centage  difference  between  the  elastic  and  crushable  foam  models  is  as  high 
as  701.  The  basic  implication  from  this  increased  dissipative  layer  thick¬ 
ness  is  that  (as  one  would  expect)  an  enlarged  energy  absorbing  layer  gives 
better  protection  to  the  interior  layers. 

The  difference  between  the  bulk  sound  speed  in  the  tables  and  diploe 
is  very  slight  (1.6  to  1.467)  in  the  basic  formulation.  The  acoustic  impe¬ 
dance  (density  multiplied  by  sound  speed)  of  a  material  layer  is  important 
in  wave  propagation  analysis.  If  there  is  a  larger  acoustic  impedance  mis¬ 
match  between  two  layers,  then  a  greater  portion  of  the  wave  energy  is  re¬ 
flected  at  the  boundary  instead  of  transmitted.  In  order  to  investigate  a 
variation  of  this  impedance  difference,  the  analytic  computer  code  was  pro¬ 
gramed  for  a  sound  speed  ratio  of  1.6  to  1.0.  Figure  (36)  displays  the  in- 
temal  energy  results  for  an  elastic -plastic  diploe  with  a  large  acoustic 
impedance  mismatch.  Again,  as  expected,  the  energy  transmitted  to  the  brain 
is  somewhat  reduced.  The  parameter  stud)'’  is  not  detailed  enough  at  this 
point  to  make  a  definitive  judgment  as  to  whether  the  diploe  thickness  para¬ 
meter  or  diploe  bulk  sound  speed  parameter  is  more  effective  in  reducing 
the  enei'gy  transmitted  to  the  biain.  Figure  (1/)  displays  the  results  of 
both  enlarging  and  decreasing  the  wave  speed  in  the  diploe  region.  With  both 
modifications  included,  there  is  now  a  percentage  difference  as  large  as  11% 
between  the  elastic  and  elastic-plastic  diploe  models.  Due  to  the  slower 
wave  speed  and  larger  distance  to  traverse,  there  is  now  an  appreciable 
phase  difference  between  the  two  curves.  The  dashed  line  indicates  that 
peak  energies  are  decreasing  with  time.  A  straight  line  on  this  semi-logar¬ 
ithmic  scale  suggests  that  a  logarithmic  energy  decay  occurs. 
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The  crushable  foam  diploe  model  yielded  results  which  were  less  pre¬ 
dictable  than  those  of  the  previous  examples.  When  the  acoustic  impedance 
decreased,  the  energy  in  the  brain  increased  (see  Figure  (18)) .  As  the 
extra  energy  was  reflected  back  into  the  crushable  diploe  region,  the 
pressure  in  the  diploe  was  increased  by  25%.  While  the  material  is  crush¬ 
ing  to  its  solid  condition,  the  sound  speed  is  changing  from  the  value  for 
foam  to  the  value  for  solid  material.  An  increased  pressure  causes  the 
foam  to  rapidly  crush  towards  the  solid  sound  speed  condition.  The  over¬ 
all  significance  is  that  the  mathematical  non-linearities  of  the  hydro- 
dynamic  foam  model  suggest  that  a  sound  speed  parameter  optimization  exists 
for  foam  energy-absorbing  layers.  For  the  elastic -plastic  diploe  model, 
the  total  plastic  work  done  in  the  diploe  can  be  computed .  The  value  of 
plastic  work  as  a  function  of  time  is  plotted  in  Figure  (19).  The  signific¬ 
ance  of  plastic  work  is  that  it  represents  stored  energy  which  cannot  be 
transmitted  into  the  brain  as  injurious  energy. 

No  previous  one -dimensional  models  are  formulated  in  an  analogous 
manner  to  the  present  analysis.  The  lumped -parameter  models  cannot  predict 
local  stress  levels  and  the  rigid  container  model  of  Hayashi  (1969)  is 
accurate  for  soft  impacts  before  the  container  rebounds.  The  "exact"  ex¬ 
pansion  solution  of  Liu  (1970)  predicts  wave  propagation  effects,  but  the 
rigid  container  aspect  results  in  zero  pressure  at  the  center.  The  present 
wave  propagation  emanating  from  impact  cm  one  side  of  the  head  and  traveling 
back  and  forth  through  the  brain  seems  to  be  an  realistic  model  for  the 
actual  head  impact  situation. 

In  searching  for  the  highest  level  of  tensile  stress  near  the  blow  and 
opposite  to  the  blow,  there  is  a  chance  of  missing  the  true  peak  value  with 
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the  reduced  number  of  outputs.  However,  the  indication  of  trends  is  more 
significant,  at  present,  than  the  extra  information  gained  from  costly 
excess  computer  operation.  The  output  v/as  scanned  to  locate  the  tensile 
stress  peaks  near  the  site  of  the  blow  and  near  the  opposite  side.  In 
addition  to  determining  the  value  and  time  of  occurrence  for  the  first 

tensile  stress  peaks,  similar  information  was  obtained  for  the  second  »* 

<! 

peaks.  Table  4  displays  the  first  and  second  peak  values  and  their  re-  * 

spective  times  for  several  diploe  models.  All  values  listed  are  in  non- 

j 

dimensional  units  (one  unit  is  approximately  SOu  sec) .  In  a  hydrodynamic 

model,  normal  stress  is  equivalent  to  pressure.  Tensile  pressure  is  con-  , 

sidered  positive  in  Table  4. 


MAGNITUDE  AND  OCCURRENCE  OF  PEAK  TENSILE  STRESSES 


SKULL 

LAYER 

MODELS 

. . -  ~  -  - .  -  ■■  ■  —  "  - 

1st  peak  1st  peak  2nd  peak  2nd  peak 

near  opposite  near  opposite 

blow  time  blow  time  blow  time  blow  time 

xl0'2  xlO'2  xlO'2  *’  xlO’2 

.738 

n 

.662 

H 

.668 

7.8 

.432 

ELASTIC- 

PLASTIC 

.566 

KB 

.601 

3.5 

.578 

s 

.373 

8.2 

CRUSHABLE 

.526 

3.0 

.624 

3,6 

.626 

m 

.432 

8.3 

ELASTIC: 

PIASTIC 

.588 

g 

.594 

D 

.617 

8.6 

.385 

9.5 

*With  enlarged  diploe  and  large  impedance  mismatch. 

TABLE  4 

PEAK  TENSILE  STRESSES  NEAR  THE  SITE  OF  THE 
BLOC  AND  OPPOSITE  THE  SITE  OF  'DIE  BLOW 


(Stress  and  time  are  given  in  non-dimensional  units.) 


70 


The  presence  of  relatively  high  tensile  stresses  near  the  site  of  the 
blow,  near  the  opposite  side  from  the  blow,  and  even  in  the  central  regions 
suggest  potential  sites  of  coup,  contrecoup,  and  intermediate  coup.  The 
important  fact  to  note  is  that  these  tensile  stresses  can  appear  in  the 
absence  of  a  closed  container  model.  This  fact  agrees  with  the  concept  of 
Goldsmith  (1966)  that  rarefaction  waves  reflecting  from  the  interior  skull 
surface  cause  high  tensile  stress  levels. 

The  information  shown  in  Table  4  was  obtained  from  the  second  spatial 
zone  and  next  to  last  zone  for  the  brain  layer  of  the  one -dimensional  model. 
The  sharp  rise  of  a  stress  wave  could  cause  an  absolute  peak  value  to  be 
unobserved  by  the  data  output  system.  However,  the  peak  tensile  stresses 
in  the  brain  seem  to  be  reduced  by  the  energy  absorbing  diploe  layer  models. 
There  is  also  a  general  trend  for  reduced  pressure  at  other  brain  locations. 
The  stress  levels  near  the  site  of  the  blow  appear  to  remain  relatively 
high  at  the  second  peak,  while  the  opposite  side  of  the  brain  experiences  a 
reduced  stress  level  at  the  second  peak. 

An  important  concept  in  the  potential  mechanisms  for  head  injury  is  the 
duration  of  the  injurious  force.  Liu,  et  al  (1971)  define  a  Cumulative 
Damage  Criterion  which  is  based  on  the  pressure-time  area  beyond  a  cut-off 
tensile  stress  level.  Therefore,  a  possible  design  criterion  for  protective 
head  devices  is  not.  only  reduction  of  energy  levels  and  peak  tensile  stress, 
but  also  control  of  the  duration  of  the  peak  stresses.  The  present  computer 
output  is  not  adequate  to  accurately  analyze  the  peak  duration  times.  A 
scanning  of  the  data  yields  the  approximate  durations  listed  in  Table  5.  It 
is  to  be  noted  that  a  peak  cut-off  level  of  2.0  x  10  ^  non-dimensional  units 
of  stress  was  used  to  determine  the  values  for  the  duration  times.  The 
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practical  result  of  using  the  above  cut-off  level  was  that  the  entire  ten¬ 
sile  stress  wave  was  observed.  Table  5  indicates  only  duration  and  not 
tensile  stress  level;  therefore,  the  results  are  not  directly  analagous  to 
the  Cumulative  Damage  Criteria.  The  site  of  potential  contre-coup  injury 
is  exposed  to  tensile  stress  for  shorter  periods  than  the  potertial  coup 
injury  location.  The  slightly  longer  durations  with  the  crushable  and 
elastic-plastic  models  indicate  the  general  smoothing  and  spreading  of  the 
wave  which  is  the  result  of  including  energy  absorbing  mechanisms. 


TABLE  5  -  Peak  Tensile  Stress  Duration 


DIMENSIONLESS  TIME  DURATION 


SKULL 

LAYER 

4DDELS 

— 

— 

1st  peak 
near  blow 

1st  peak 
opposite 
blow 

2nd  peak 
near  blew 

2nd  peak 
opposite 
blow 

elastic 

2.3 

.7 

1.4 

1.1 

elastic- 

plastic 

2.0 

.8 

1.6 

.8 

crushable 

2.2 

.9 

1.7 

.8 

clast ic- 
slastic 

2.5 

1.0 

2.0 

1.3 

*Wjth  enlarged  diploe  and  large  impedance  mismatch. 
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Hie  results  of  the  rate  dependent  brain  model,  for  all  diploe  models, 
are  shown  in  Figure  20.  The  shape  of  the  curves  is  similar  to  the  shape 
of  the  rate  independent  brain  model  curves.  The  peak  internal  energy 
values  for  all  diploe  models  are  approximately  equivalent  for-  the  two 
brain  models.  The  lower  portion  of  the  internal  energy  curves  for  the 
energy  absorbing  diploe  models,  with  the  viscous  brain  layer,  was  slightly 
higher  than  that  of  the  inviscid  brain  model.  There  are  no  conclusive 
explanations  as  yet  for  this  increased  energy.  A  potentially  larger 
acoustic  impedance  mismatch  could  lead  to  extra  energy  storage  in  the 
brain  region.  Because  the  viscous  model  is  first  order  accurate  (compared 
to  second  order  accuracy  for  all  other  constitutive  models) ,  an  increased 
number  of  time  cycles  is  required  to  attain  meaningful  results.  Despite 
the  attanpt  to  achieve  a  uniform  degree  of  accuracy,  the  rate-dependent 
model  yie.ds  energy  balance  errors  as  large  as  101.  This  error  could 
have  produced  the  higher  energy  found  in  the  viscous  model.  In  general, 
the  viscous  brain  model  does  not  seem  to  differ  greatly  from  the  inviscid 
model,  but  further  analysis  is  required  to  more  accurately  describe  the 
rate-dependent  brain  effects. 


TWo  Dimensional  Analysis 

Four  spherical  model  analyses  were  solved  as  follows:  (1)  elastic 
single  layered  skull  and  hydrodynamic  brain;  (2)  elastic  single  layered 
skull  and  elastic  (shear)  brain;  (IS)  layered  skull  with  elastic-plastic 
diploe  and  hydrodynamic  brain,  and  (4)  layered  skull  with  crushaMe  diploe 
and  hydrodynamic  brain.  While  some  skull  data  is  presented,  the  primary 
effort  is  concentrated  on  determining  potentially  dangerous  mechanical 
stress  conditions  that  exist  in  lie  brain.  The  first  and  h.ai  model 
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analyses  were  solved  for  8.0  units  of  dimensionless  time,  while  the  second 
and  third  are  only  solved  for  3.5  units. 

The  large  volume  of  data  generated  by  the  "TOODY-rV”  computer  code  can 
only  be  interpreted  by  considering  selected  graphic  displays.  Three  forms 
of  computer  drawn  plots  are  utilized.  Figure  21  depicts  a  polar  coordinate 
system  for  identifying  the  locations  at  which  pressure  values  are  plotted. 

The  0°  line  actually  represents  the  average  stress  values  in  the  zone  adjacent 
to  the  symmetry  line.  The  180°  line  is  defined  in  the  same  way.  By  consider¬ 
ing  a  fixed  angle,  the  pressure-time  curve  for  a  given  radius  may  be  drawn. 

It  should  be  noted  that  in  all  plots  tensile  pressure  is  considered  positive. 

To  present  an  overall  perspective  of  the  pressure  at  several  radii  along  a 
given  angle,  a  bidder.,  line  plot  is  used  which  shows  only  the  maximum  tensile 
pressure  that  is  visible  when  six  radii  are  drawn  with  an  effective  three- 
dimensional  depth  into  the  paper.  A  final  type  of  graph  displays  pressure 
ever  tire  entire  brain  region  a':  a  given  time.  Each  hemisphere  (the  symmetric 
image  is  equivalent)  contain.?  144  values  of  relative  pressure  from  +9  to  -9. 
Table  6  desert..  ;s  the  cor.  "ersion  from  dimensionless  pressure  to  the  relative 
scaie  used  in  +he  p.t.-ssure  profiles. 

The  singi*.'  jaycred  elastic  skull  and  hydrodynamic  brain  model  is  the 
basic  model  that  serves  as  a  reference  for  other  models  in  this  paper  and  as 
a  comparison  with  previous  analyses.  Even  this  simple  model  is  very  complex 
in  terms  of  physically  c>p  aining  the  wave  propagation  results.  All  pressure 
and  stress  values  are  tie  composite  result  of  several  mechanical  phenomena. 

As  the  primary  wave  passes  through  the  skull  anJ  cult ”b  the  brain,  a 
stress  wave  travels  around  the  skull  becoming  a  new  point  source  for  generating 
waves  inward  towards  the  brain.  Because  the  jmjiactor  only  strikes  a  small 
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portion  of  the  skull,  the  regions  adjacent  to  the  indentation  bulge  outward 
giving  rise  to  tensile  relief  waves.  Subsequent  to  the  initial  in-bending 
of  the  shell  near  the  symmetry  axis,  the  shell  snaps  back  to  produce  rare¬ 
faction.  waves  from  outbending.  When  any  of  these  waves  travelling  in  one 
medium  impinges  on  a  boundary  region  of  another  medium,  then  all  the  stress- 
optical  laws  of  reflection  and  refraction  are  obeyed.  The  most  significant 
of  these  wave  phenomena  is  the  change  from  compression  to  rarefaction  that 
occurs  when  a  wave  moving  in  a  medium  of  high  acoustic  impedance  reflects 
from  the  surface  of  a  lower  acoustic  impedance  material.  Also  the  skull 
tends  to  accelerate  its  mass  center  before  the  brain  starts  to  move;  there¬ 
fore,  on  the  side  opposite  from  the  impactor  there  is  an  additional  source 
of  tensile  force.  Because  the  primary  objective  of  this  analysis  is  to 
evaluate  damage  potential  in  the  Train,  there  is  only  a  secondary  effort 
to  ascertain  the  source  of  the  resultant  wave  propagation. 


Impact  to  the  elastic,  skull  and  hydrodynamic  brain  model  produced 
regions  of  high  tensile  pressure  in  the  brain  near  the  blow,  on  the  opposite 
side  fron  the  blow,  and  at  several  locations  in  between.  Using  the  approach 
that  brain  tissue  fails  by  excessive  tensile  pressure  (see  Benedict  (1970) 
and  Iingin  (1969)),  these  regions  of  rarefaction  could  represent  the  mechanism 
for  coup,  centre -coup,  and  intermediate  coup,  respectively.  TV/o  of  the 
three-dimensional  perspective  graphs  show  the  pressure  results  along  the 
symmetry  axis,  figures  22  and  23  display  the  pressure  -  time  curves  for  six 
ladti  along  the  0°  and  130°  lines  in  the  brain.  ,;or  times  up  to  3,!j 
dimensionless  units,  1)/:  largest  tensile  pressure  appears  at  the  outer  mdius 
along  the  IBU’  line.  Lesser  rarefaction  peaks  occur  at  several  points  through¬ 
out  the  brain.  At  angles  further  away  from  the  symmetry  line,  the  picssure-  is 


skull  two -dimensional  model. 
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generally  weaker.  Figure  24  shows  the  reduced  pressure  levels  that  exist 
along  the  96°  line. 

The  exact  values  of  pressure  are  readily  observed  in  the  standard 
pressure  curve  giver  in  Figure  25.  Because  physiological  head  injury. cannot 
be  attributed  to  a  definite  mechanism,  there  are  several  possible  pressure 
injury  phenomena  that  can  be  observed  from  Figure  25.  Two  types  of  peak 
tensile  pressures  appear  at  the  contre- coup  site.  Sharp  spike  peaks  occur 
at  times  of  approximately  2.0  and  2.8  units  of  time,  and  a  broader,  but 
lower  level  peak  occurs  at  about  4.0  units  of  time.  If  peak  level  (regard¬ 
less  of  duration]  causes  tissue  and  vessel  damage,  then  the  earlier  peaks 
could  cause  injury.  Also  the  rapid  fluctuation  between  highly  positive  and 
negative  pressure  would  tend  to  enhance  a  mechanism  that  depends  on  instan¬ 
taneous  failure  at  peak  levels.,  Liu,  et  al  (1971)  proposed  that  the  time 
averaged  pressure  (analogous  to  impulse)  was  the  critical  injury  mechanism. 
Using  the  entire  area  of  the  tensile  pressure  pulse  at  2.8  and  4.0  dimension¬ 
less  times,  the  broader  pressure  peak  (at  die  latter  time)  has  a  greater  area. 
Figure  26  shows  the  entire  pressure  field  .in  the  brain  at  a  time  of  2.8  units. 
The  wedge  shaped  region  (only  half  is  shown  in  the  symmetric  hemisphere)  of 
high  tensile  pressure  opposite  fron  the  blew  and  near  the  suds  line  looks 
very  similar  to  th-  typical  pathological  finding  in  head  injury  victims  (c.g. 
Lindenberg  and  1  rcytag  (1957)).  Wedge  shaped  regions  were  also  observed  with 
the  broader  pulse. 

The  intermediate  location  (dimensionless  radius  of  .17)  is  subjected  to 
several  large,  spike  shaped  tensile  pulses.  'Ihc  peak  values  arc  not  as  high 
as  at  the  outer  radius,  and  the  pressure  profile  shown  in  Figure  27  indicates 
tliut.  this  tensile  pressure  concentration  is  loss  diffuse  than  at  the  outer 
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radius.  As  indicated  in  Figure  22,  the  rarefaction  peaks  near  the  blow  ; 

J* 

are  much  veakcr  and  much  less  diffuse  than  those  opposite  the  blow.  i 

3 

The  elastic  skull  results  described  above  compare  favorably  with  the  ( 

% 

analytic  solutions  of  similar  problems  with  impulsive  loads  as  solved  by  \ 

a 

Engin  (1969),  Iiu,  et  al  (1971)  and  Chan  (1972).  The  time  of  occurrence 

<C 

and  general  levels  of  the  peak  tensile  stresses  for  contre-coup  and  inter-  t 

mediate  coup  are  in  close  agreement  with  the  results  of  these  authors.  The 

rapid  transition  from  compressive  to  tensile  pressures  is  discussed  by  Liu, 

et  al  (1971).  This  latter  article  also  discussed  the  large  amplifications  ■ 
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of  input  pressures  that  result  in  high  peak  pressures  (1.8  x  10  dynes/an 
in  the  present  model)  when  dissipating  mechanisms  are  not  included. 

The  values  of  pressure  and  shear  stress  in  the  skull  provide  some  in 

sight  into  possible  regions  of  skull  fracture.  Evans  (1957)  describes  skull 

fracture  as  being  directly  related  to  tensile  stresses  at  and  near  the  impact 

site.  Figure  28  displays  the  ski'll  pressure  (average  of  the  normal  stresses)  ; 

for  an  outer,  middle,  and  inner  radius  of  the  skull  along  the  0°  line.  The  ; 

/ 

outer  radius  shows  the  high  compressive  pressure  that  is  the  result  of  decreased  ; 

curvature  from  inbending  plus  strong  compression  waves  from  the  impact or . 

The  inner  radius  would  be  in  tension  because  of  increased  curvature,  but  the 
impact  compression  waves  produce  a  combined  mildly  compressive  initial  pressure. 

At  later  times  the  snap  back  yields  high  rarefaction  pressures  in  the  outer 
shell.  High  tensile  pressures  in  the  outer  skull  appear  earlier  along  the 
24°  line.  Shear  levels  are  also  high  near  the  blow  and  close  to  the  symmetry 
axis  (see  Figure  29).  The  general  result  of  potential  sites  of  skull  failure 
existing  in  an  arc  from  0°  to  approximately  35°  are  in  -g* cement  with  the 
discussions  of  Fngin  (1969).  Except  for  some  small  pressure  fluctuations  at 
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the  side  opposite  the  impactor,  the  skull  stresses  in  the  intermediate 
angles  are  relatively  weak.  The  regions  of  high  tensile  pressure  in  the 
FkL-'.i  could  possibly  be  related  to  the  Maximum  Strain  Index  (McElhaney, 
et  al  (1971))  in  the  skull  when  the  problem  is  extended  to  longer  solu¬ 
tions  . 

The  oscillation  between  inbending  and  outbending  implies  an  apparent 
shell  frequency  of  approximately  2000  Hz.  This  value  is  over  twice  the 
experimental  value  of  Stalnaker,  et  al  (1970)  and  Hodgson  and  Patrick  (1963). 
The  explanations  for  this  discrepancy  are  that  the  theoretical  results  are 
based  on  dynamic  transient  effects  and  that  the  experimental  head  mass  in¬ 
cludes  the  ext^a  mass  of  the  facial  bones. 

The  elastic  skull  and  elastic  brain  model  results  provide  an  evaluation 

of  the  shear  stress  that  may  be  induced  from  an  axisyirmetric  translational 
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blow.  The  assumed  dynamic  shear  modulus  of  10  dynes/cm  is  less  than  one 
tenth  the  bulk  modulus;  therefore,  the  resultant  shear  stress  values  are 
less  than  one  tenth  of  the  pressure  values.  For  the  elastic  brain  model, 
the  pressuie  values  are  computed  as  the  average  of  the  normal  stresses. 

Figure  30  shows  that  the  general  form  of  the  pressure-time  curves  is  similar 
to  the  results  for  the  hydrodynamic  model.  The  peak  pressures  are  slightly 
reduced  because  some  of  the  energy  is  incorporated  into  the  shear  stress. 

Shear  values  for  the  x-z  plane  (that  are  available  from  the  "TOODY-IV” 
code)  are  shown  in  Figure  31.  The  168°  line  is  used  because  these  values 
are  higher  than  those  at  the  180°  line  (note  that  the  scale  used  in  one- tenth 
of  tne  hydrodynamic  results  scale) .  There  is  no  definite  reason  for  this 
angular  shift  in  the  peak  shear  levels.  Possibly,  the  typically  slower  speed 
cf  shear  waves  has  induced  a  phase  lag. 
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The  shear  and  pressure  profiles  at  the  dimensionless  time  of  2.77  are 
shown  in  Figures  32  and  33.  The  high  simultaneous  values  of  shear  and  tensile 
pressure  suggest  chat  the  combined  effect  of  these  two  mechanical  forces 
might  produce  tissue  or  vessel  damage.  Typical  shear-normal  failure  plots 
for  most  materials  show  that  adding  shear  enhances  breakage  at  lower  normal 
stress  levels.  The  use  of  a  high  dynamic  shear  modulus  (as  in  the  case  of 
dilute  polymer  gels)  also  implies  a  brittle  "glassy"  state  exists  that  would 
result  in  a  fracture  mechanism  injury.  In  light  of  the  experimental  and 
pathological  evidence  of  shear  produced  damage,  the  shear-tensile  pressure 
injury  concept  seems  to  provide  a  reasonable  and  unifying  explanation  of 
tissue  and  vessel  failure. 

The  existence  of  a  shear  and  normal  stress  at  a  given  point  means  that 
there  is  a  principle  plane  (at  that  point)  on  which  the  normal  stress  is 
maximized.  Near  the  outer  radius  along  the  168°  line  at  the  2.77  units  of 
time,  the  stress  conditions  are  T**  =  .79  x  10~3,  Tzz"=  .83  x  10’3,  and  T*2  = 
.114  x  10  3  dimensionless  units  of  stress.  A  plane  exists  on  which  the  normal 
stress  if  .93  x  10  J,  At  other  locations  with  lower  shear  levels,  the  in¬ 
crease  in  principle  stress  is  not  as  large. 

The  multi-layered  skull  and  hydrodynamic  brain  models  determine  the 
effects  of  a  central  skull  region  of  lesser  acoustic  impedance  and  of  energy 
absorption  in  that  region.  Because  it  was  desired  to  accentuate  the  energy 
absorbing  mechanisms  (the  same  goal  can  be  accomplished  with  a  higher  energy 
impact),  the  value  of  the  yield  stress  and  initial  crushing  pressure  are  one- 
tanth  of  the  static  experimental  values  of  Melvin,  et  al  (1970).  Skull 
pressures  in  all  of  these  analyses  were  higher  than  with  an  elastic  skull. 
Figure  34  shows  the  skull  pressure  adjacent  to  the  impactor  along  the  symmetry 
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axis  for  a  crushable  diploe  model. 

The  elastic-plastic  diploe  model  limits  the  deviators  as  described  by 
equation  39.  Figure  35  displays  the  generally  reduced  tensile  pressures  in 
the  brain  that  result  from  impact  to  this  model.  The  peak  contre-coup  rare¬ 
faction  pressure  is  reduced  by  one  half  compared  to  the  elastic  model.  There 
is  still  substantial  rapid  positive  to  negative  pressure  fluctuation  that 
indicates  the  pressure  waves  are  not  greatly  smoothed.  If  a  more  efficient 
yield  stress  mechanism  (such  as  work  or  strain  hardening)  was  experimentally 
determined  to  exist  in  the  diploe,  then  even  greater  pressure  reduction  and 
more  smoothing  could  be  obtained. 

The  crushable  foam  diploe  model  is  a  more  efficient  energy  absorber  than 
the  elastic-plastic  model,  and  it  is  specifically  designed  to  reduce  the  peak 
levels  of  pressure  pulses.  Therefore,  the  protection  afforded  by  this  model 
yielded  the  most  significant  reduction  of  potentially  dangerous  mechanical 
forces . 

To  provide  some  continuity  with  previous  experimental  investigations,  a 
crushable  diploe  model  was  analyzed  utilizing  the  initial  crushing  pressure 
and  final  crushnp  pressure  as  evaluated  by  Melvin,  et  al  (1970)  and  McElhaney, 
et  al  (1970).  The  general  levels  and  peak  tensile  pressures  were  slightly 
lower  than  for  the  elastic-plastic  model  with  a  lower  yield  point.  No  plots 
were  constructed  for  this  model.  The  reduced  pressure  was  observed  from  the 
data  printed  uuring  the  computer  analysis. 

Figure  36  shows  the  pressure  along  the  180°  line  for  the  crushable  model 
with  a  lower  initial  crushing  value.  The  general  smoothing  of  the  curves 
represents  a  drastic  modification  compared  to  the  elastic  model  results.  The 
highest  rarefaction  occurs  at  about  2.9  units  of  time.  With  the  elastic  skull 
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model  there  was  a  peak  at  the  time  of  2.0,  but  a  larger  value  appeared  at  2.8 
units.  The  crushable  analysis  peak  was  less  than  one  fourth  of  that  for  the 
elastic  analysis.  The  pressure  profile  at  the  time  of  peak  rarefaction  (see 
Figure  37)  shows  a  wedge  shaped,  very  diffuse,  and  lower  valued  concentration 
of  tensile  pressure  at  the  contre-coup  site.  The  combination  of  acoustic  im¬ 
pedance  mismatch  plus  strong  energy  absorption  in  the  diploe  region  result  in 
substantial  protection  against  head  impact  injury. 

The  internal  energy  in  the  brain  region  was  not  available  as  magnetic 
tape  output  from  the  computer  analysis.  Therefore,  the  hand  plotted  internal 
energy  versus  time  curve  appears  in  Figure  38.  The  general  periodic  peaking 
was  not  as  regular  as  for  the  one -dimensional  analysis.  While  stress  waves 
carry  energy  into  and  out  of  the  brain  region,  there  are  specific  times  at 
which  the  spherical  shape  produces  strong  reinforcement  peaks.  The  higher 
initial  crushing  value  model  transmitted  less  energy  to  the  brain  than  the 
elastic -plastic  model  with  a  reduced  yield  stress.  The  lower  crushing  value 
foam  model  yielded  the  least  energy  in  the  brain.  While  internal  energy  re¬ 
presents  positive  as  well  as  negative  stress,  it  is  interesting  to  note  that 
the  elastic  model  shows  peaks  at  the  dimensionless  times  of  2.0  and  2,8  and 
the  lower  crushing  value  foam  model  shows  no  peak  at  2.8.  The  pressure  plot 
observations  of  rarefaction  peaks  in  the  brain  with  the  addition  of  more 
energy  absorption  in  the  diploe  are  closely  matched  by  the  reduced  internal 
energy  results.  This  indicates  the  general  usefulness  of  this  method  for 
comparing  protection  from  different  head  models. 
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Elastic  Diploe 

Crus! able  Diploe  (low  crushing  point) 
Crushable  Diploe'  (high  crushing  point) 


Figure  38.  Internal  energy  in  brain  of  two-dimensional 
analyses . 
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Experimental  Analysis 

Several  selected  photographic  results  from  the  dynamic  photo-elastic  ex¬ 
periments  are  shewn  in  Figure  39  to  Figure  45.  The  dark  fringe  wave  front 
is  traveling  from  right  to  left  in  all  frames,  and  the  time  difference  be¬ 
tween  frames  is  indicated  under  the  photographs  (if  applicable).  The  dark 
circle  at  the  boundary  between  inner  disk  and  outer  ring  and  around  a  few 
diploe  holes  was  the  result  of  a  pre- stress  from  the  machining  and  casting 
processes.  Despite  this  flaw  in  the  model  making  techniques,  the  stress 
wave  propagation  is  clearly  observable.  The  original  photographs  are  Polaroid 
prints  with  each  frame  being  about  .80  cm  by  1.50  an.  In  certain  frames  with 
poor  lighting,  the  blown  up  figure  is  unfortunately  not  a  high  quality  print. 


Plane  waves  entering  the  small  holes  of  the  table  models  remained  almost 
plane,  and  completely  continuous  while  traversing  the  holes  and  also  emerged 
as  a  plane  wave  (as  shown  in  Figures  39  and  40) .  The  pulse  size  of  the  main 
pressure  front  can  be  approximated  by  multiplying  the  pulse  duration  (1 Op sec, 
as  determined  from  preliminary  studies)  times  the  wave  speed  in  PSM-5  (1830 
cm/sec).  The  resulting  pulse  width  of  1.83  cm  is  much  larger  than  each  indi¬ 
vidual  hole  depth  so  that  the  wave  can  start  to  re -group  before  the  next 
line  of  holes  is  encountered  (one  diameter).  Previous  studies  (Rose,  Mortimer, 
arid  Chou  (1973))  indicated  that  for  pulse  widths  at  least  as  large  as  the 
hole  depth  the  wave  front  can  regroup  in  two  characteristic  lengths.  Character¬ 
istic  length  is  defined  as  the  hole  dimension  parallel  to  the  wave  front.  It 
should  be  noted  that  the  lining  up  of  small  holes  parallel  to  the  wave  front 
helps  to  maintain  a  wave  with  slightly  rippled  form,  but  does  little  to  aid 
in  obtaining  a  continuous  wave.  The  small  characteristic  length  of  each  hole 
provides  the  primary  enhancement  for  re-grouping  as  a  continuous  wave.  Also 


i 


l 


i 


102 


Figure  39.  Photoelastis  stress  waves  in  table  region  with 
PL- 2  filler. 


exiting 


entering 


Figure  40.  Photoeiastic  stress  waves  entering  and  exiting  table 
region  with  Solithane  113  filler. 
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the  entire  table  region  is  thinner  in  depth  than  the  pulse  width;  therefore, 
this  total  section  acts  as  an  effective  inclusion  vjhich  can  be  traversed  as 
a  continuous  and  plane  wave. 

PL-2  and  Solithane  113  models  yielded  approximate  wave  speeds  in  the 
small  hole  region  of  1710  and  1690  cn/sec  respectively.  These  values  are 
nearly  equal  (despite  the  difference  in  filler  properties)  and  are  closer 
to  the  PSM-S  wave  speed  than  to  the  filler  wave  velocity.  The  similarity 
in  wave  speed  and  wave  form  results  of  the  two  table  models  suggests  that 
the  geometrical  considerations  of  small  individual  hole  size,  closely  packed 
arrangement  and  thin  overall  depth  are  more  significant  than  the  filler 
material  utilized. 

The  larger  holes  of  the  diploe  region  had  a  more  pronounced  effect  on 
the  wave  motion.  With  the  PL- 2  filler  the  wave  front  was  highly  irregular, 
but  continuous.  However,  the  Solithane  113  model  produced  a  non-continuous 
wave  that  energed  as  distinct  rays  emanating  from  the  solid  portions  of  the 
model.  The  larger  characteristic  length  of  each  hole  did  not  enhance  re¬ 
grouping  before  the  next  cavity  was  met.  Also  the  entire  diploe  region  was 
larger  than  the  input  pulse  width.  The  general  effect  was  that  only  in  the 
case  of  a  small  acoustic  impedance  mismatch  could  the  wave  remain  continuous. 
These  wave  propagation  results  can  be  observed  in  Figures  41  and  42.  The 
concentration  of  stress  in  the  solid  portions  of  the  diploe  indicate  a  region 
of  early  failure  as  suggested  by  Melvin,  et  al  (1970). 

Because  of  the  large  depth  of  the  entire  diploe  section  it  was  impossible 
to  observe  a  wave  entering  and  exiting  the  diploe  during  the  I2ysec  period  of 
a  single  test.  By  noting  the  time  duration  between  a  wave  exiting  the  large 
holes  (in  the  Solithane  113  model)  and  a  wave  passing  through  a  similar 
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4  psec 


Figure  41.  Photoelastic  stress  waves  in  diploe  region  with  PL-2 
filler. 


ngure  42.  Photoelastic  stress  waves  in  diploe  region  with  Solitham 
11a  filler. 
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location  on  the  reference  plate  (opposite  side  of  the  exploding  wire),  the 
overall  average  speed  in  both  regions  of  the  skull  model  is  computed  as  1S25 
an/sec.  This  overall  value  suggests  that  the  speed  in  the  diploe  is  slower 
than  in  the  table  region. 

Ail  of  the  exploding  wire  tests  are  in  the  elastic  range  because  there 
are  no  permanent  deformations  in  the  table  or  diploe  regions.  Therefore, 
the  experimental  comparisons  with  the  theoretical  analysis  are  applicable 
for  the  initial  time  period  before  crushing  can  begin.  The  slower  wave 
speeds  in  the  diploe  region  are  physically  predictable  and  in  agreement  with 
the  findings  of  the  West  Virginia  Biomechanics  Laboratory  (1970) .  These  re¬ 
sults  add  confidence  to  the  relative  wave  speed  values  employed  in  the 
theoretical  analysis. 

The  exact  nature  of  wave  propagation  is  of  course,  strongly  dependent 
on  the  relationship  between  hole  size  and  pulse  size.  While  the  plastic 
skull  model  is  proportional  to  the  human  head,  the  hole  size  to  pulse  size 
ratio  for  the  model  is  higher  than  for  direct  impact  to  ai  actual  head.  The 
photoelastic  results,  using  PL-2  which  is  closer  to  the  fluid  filled  diploe 
properties,  can  be  extended  to  yield  a  qualitative  approximation  of  the 
wave  form  during  head  impact.  The  microscopic  cavities  of  the  table  regions 
probably  have  little  effect  on  the  wave  motion  other  than  to  slightly  reduce 
the  wave  speed.  The  larger  spaces  in  the  diploe  (but  smaller  than  the 
plastic  model  holes)  yield  a  continuous  but  not  perfectly  planar  wave  foim 
with  a  lower  wave  velocity  than  the  table  region.  Because  the  non-planar 
motion  in  the  diploe  will  reform  into  a  plane  wave  after  traveling  through 
the  inner  table,  the  ..mall  ripples  in  the  wave  form  are  considered  secondary 
compared  to  the  crushing  effects.  Additionally,  it  would  be  theoretically 
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impossible  to  incorporate  all  of  the  inhomogeneous  diploe  inclusions  in  the 
analysis . 

The  disk  and  concentric  ring  photoelastic  head  model  provided  direct 
observation  of  stress  waves  traveling  through  the  brain  and  around  the  skull 
model  analogs.  A  mathematical  analysis  of  this  plastic  head  model  (with  PL-2 
disk)  was  performed  using  the  "TOCDY-TV”  computer  code.  A  cylindrical 
theoretical  geometry  was  used  (similar  to  the  shape  shown  in  Figure  4-a) 
that  does  not  allow  for  possible  edge  effects  present  in  a  thin  plastic 
model.  Except  for  this  infinite  thickness,  the  cross-sectional  configura¬ 
tion  is  the  same  for  both  models.  The  early  time  results  for  each  analysis 
(least  affected  by  edge  effects)  have  been  compared,  and  the  correlations 
are  in  close  agreement . 

The  wave  form  for  both  analyses  had  similar  characteristic  shapes.  As 
the  wave  front  started  to  enter  the  central  disk,  the  wave  had  traveled 
approximately  45°  in  the  outer  ring  of  each  model.  The  curved  wave  form  in 
the  shell  oi  both  models  indicated  the  slower  wave  velocity  due  to  contact 
with  lower  acoustic  impedance  materials  at  the  inner  and  outer  shell  radii. 

A  slight  advancement  of  the  wave  front  along  the  symmetry  axis  as  compared 
to  other  locations  showed  the  geometric  dispersion  that  exists  in  the  disk 
region.  As  pressure  starts  to  be  generated  inward  from  the  shell  (as  the 
wave  passes  around  the  cuter  ring) ,  the  wave  front  in  the  disk  attains  a 
crescent  shape  (see  Figure  43).  When  the  wave  reaches  the  opposite  side 
from  the  exploding  wire,  there  is  a  strong  concentration  of  stress  along 
the  outer  radius  of  the  disk.  Figure  44  shows  the  distinctive  horseshoe 
shape  that  resulted  with  the  Solithane  113  disk.  The  acoustic  impedance 
mismatch  is  so  large  that  the  wave  is  almost  complete!/  around  the  cuter 
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Figure  43.  Photoelastic  stress  waves  in  head  model  with  PL-2 
brain  region. 


Figure  44.  Photoelastic  stress  waves  in  head  model  with  Solithane 
113  brain  region. 
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ring  before  the  wave  enters  the  disk;  therefore,  pressure  generated  inward 
from  the  ring  has  gone  through  a  larger  angle  than  with  the  PL- 2  disk. 

The  dark  stress  fringes  are  representative  of  an  isochromatic  line 
which  implies  a  constant  difference  between  two  principle  stresses.  The 
theoretical  results  are  displayed  as  a  plot  of  the  principle  stress  differ¬ 
ences  versus  position  at  4y  sec  intervals.  In  order  to  shew  the  principle 
stress  difference  profiles,  relative  values  indicated  are  based  on  a  linear 
scale  from  1  to  9.  The  actual  leading  and  trailing  edges  of  the  wavefront 
include  several  values  lower  than  1  on  the  relative  scale.  These  lower 
level  principle  stress  differences  are  all  listed  as  0. 

Figure  45  shows  the  photographic  experimental  and  calculated  theoretic¬ 
al  wave  forms  for  two  similar  times.  The  mathematical  wave  form  is  sketched 
by  using  the  two  highest  values  in  the  outer  ring  and  the  highest  value  in 
the  disk.  Lower  isochrcmatic  theoretical  lines  are  found  in  the  disk  region 
because  of  the  lower  shear  modulus.  Both  symmetric  halves  of  the  theoretical 
results  in  the  ring  and  disk  (not  the  short  support  stmt)  are  presented. 

The  primary  wave  fronts  in  both  mod  ?ls  have  a  similar  form.  The  probable 
causes  of  discrepancies  between  the  two  models  are  imperfections  *n  the 
molded  and  machined  plastic  model,  inaccurate  mathematical  modeling  of  the 
pressure  pulse,  and  incorrect  material  property  values  for  the  theoretical 
analysis.  The  general  agreement  between  the  models  provides  qualitative 
confidence  in  the  theoretical  solution  technique. 


Figure  45.  Correlation  of  experimental  and  theoretical  ring  and  con¬ 
centric  disk  head  model. 
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VI  -  SUMMARY  AND  CONCLUSIONS 

There  is  a  lack  of  high  frequency  experimental  data  to  describe  the 
mechanical  properties  of  the  head.  The  present  analysis  is  based  on  the 
best  available  test  data  and  several  approximations  and  modifications  re¬ 
quired  to  evaluate  certain  physical  effects.  While  the  analysis  does  not 
include  every  layer  of  material  in  the  head  and  the  constitutive  proper¬ 
ties  assumed  are  not  perfect,  it  represents  a  reasonable  approach  for 
theoretically  studying  several  mechanisms  that  can  strongly  influence  head 
impact  damage. 

The  ,rWONDY-IIIa"  and  "TOODY- IV  computer  codes  have  been  developed 
and  verified  experimentally  over  a  number  of  /ears.  Internal  error  checks 
supply  a  continuous  check  of  accuracy  and  stability.  The  two-dimensional 
elastic  skull  results  bear  strong  resemblance  to  the  analytic  solutions  of 
Engin  (1969)  and  Liu,  et  al  (1971).  Dynamic  exploding  wire  studios  of  im¬ 
pact  to  a  photoelastic  head  model  correlate  qualitatively  with  the  theoretic¬ 
al  two-dimensional  analysis  of  a  similar  problem.  In  general,  there  is 
confidence  in  the  applicability  and  accuracy  of  the  finite  difference  solu¬ 
tion  techniques. 

The  one-dimensional  analysis  provides  an  economic  (in  terms  of  computer 
costs)  method  of  predicting  the  influence  of  several  material  property  and 
size  modifications  in  a  geometrically  simplified  head  impact  model.  The 
greatest  modelling  error  in  this  analysis  is  the  lack  of  a  closed  container 
that  causes  pressure  amplification  due  to  the  spherical  shape.  Various 
diploe  models  indicate  the  effects  of  impact  to  regions  with  a  thicker  diploe 
and  possibly  modified  dynamic  acoustic  properties.  A  potential  parameter 
for  optimization  of  energy  absorption  is  analysed  for  the  crushable  foam 
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diploe  model.  Rarefaction  waves  are  generated  solely  from  acoustic  mis¬ 
matches  at  material  interfaces. 

The  two-dimensional  analysis  predicts  regions  of  high  tensile  pressure 
that  correspond  to  coup,  contre-coup,  and  intermediate  coup  injuries.  Two 
types  of  pressure  pulses  generated  wedge  shaped  tension  areas  at  the  contre- 
coup  site.  Very  high  level  short  duration  spikes  could  cause  failure  by 
exceeding  a  critical  tensile  pressure  level.  A  broader  and  weaker  pressure 
pulse  appeared  that  could  cause  tissue  or  vessel  damage  by  remaining  above 
a  rarefaction  injury  level  for  a  longer  duration.  Shear  stress  was  evaluated 
in  an  elastic  brain  model.  The  combined  shear-normal  stress  levels  suggest 
that  high  tensile  stress  in  the  presence  of  a  smaller  shear  stress  is  more 
likely  to  produce  failure  than  the  shear  free  stress  conditions  in  a  hydro- 
dynamic  brain  model.  Impact  to  the  multi-layered,  skull  models  yielded  large 
reductions  in  rarefaction  pressures  in  the  brain.  The  acoustic  impedance 
mismatch  and  the  energy  absorption  (especially  in  the  crushable  foam  model) 
utilized  by  including  the  diploe  layer  of  the  skull,  provided  excellent 
protection  against  impact,  brain  damage. 


There  are  two  general  conclusions  to  be  drawn  from  the  analyses  dis¬ 
cussed  in  this  work.  First,  the  quantitative  effects  of  including  a  layer¬ 
ed  energy  absorbing  skull  and  the  possible  existence  of  a  high  dynamic 
shear  modulus  for  brain  tissue  have  significant  influence  on  the  potentially 
damaging  mechanical  forces  that  result  from  head  impact.  Second,  the 
generality  of  the  solution  techniques  readily  permits  the  extension  of  the 
analyses  to  investigate  the  importance  of  additional  layers,  new  material 
definitions,  more  complex  geometries,  or  different  types  of  impact. 
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VH  *  REOCWENDATIONS 

Three  areas  for  future  research  based  on  the  findings  of  the  present 
work  are  as  follows:  (1)  modifications  to  the  head  impact  model,  (2)  head 
protective  device  models,  and  (3)  improvements  in  the  analysis  techniques. 

The  head  inpact  modelling  could  be  extended  in  several  ways.  Improved 
experimental  high  frequency  evaluation  of  head  component  material  proper¬ 
ties  and  mechanisms  is  a  primary  area  fox  head  impact  research. 

Two  otb  .lead  layers  should  be  added  to  the  present  model.  Inpact  trans- 
mittc  i.'  ’ough  the  scalp  might  yield  a  more  gradual  input  energy  pulse. 

The  influence  of  including  the  meningeal  region  (or  a  mathematical  slippage 
at  the  inner  table-brain  boundary)  should  be  determined.  Inroactor  modifica¬ 
tions  could  include  using  a  combined  translation  and  surface  vraction  load 
(as  described  by  Qian  (1971)  or  varying  the  size  arid  .veloci  tv  of  the  impact- 
or.  A  final  change  in  the  analysis  would  be  to  study  the  effects  of  a  more 
complex  geometry  such  as  an  ellipsoidal  shell.  The  present  modelling 
techniques  can  be  readily  adjusted  to  include  the  above  considerations. 

The  crush *ble  foam  diploe  model  produced  reductions  in  potentially 
damaging  forces  in  the  brain.  This  energy  absorption  concept  can  he  extend¬ 
ed  to  include  the  addition  of  protective  helmet  layers  between  the  hnpactor 
and  outer  head  layer.  Maximum  reduction  of  energy  transmission  to  the  brain 
can  be  obtained  by  varying  material  properties  and  sizes  and  choosing  the 
proper  order  of  layering. 

Excessive  computer  costs  indicate  that  a  more  economic  analysis  would 
permit  more  parameter  modifications  and  longer  extension  of  the  problem 
times.  An  accurate  combination  of  several  layers  into  one  effective  layer 
would  simplify  the  helmsfed  head  impact  model.  Another  type  of  approxima¬ 
tion  theory  (such  as  finite  element  or  method  of  characteristics)  might  be 
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more  efficient.  Several  new  techniques  for  hybrid  solution  for  systems  of 
partial  differential  equations  can  be  considered.  A  final  approach  is  the 
use  of  mathematical  manipulation  of  shell  theories  to  reduce  the  solution 
to  a  simple  analytic  form.  The  present  finite  difference  method  will  pro¬ 
vide  a  useful  tool  for  evaluating  any  new  solution  techniques. 
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APPENDIX  A 


Complete  List  of  Finite  Difference  Equations 
For  the  Two-Dimensional  Analysis 

In  Chapter  III  -  Theoretical  Modelling,  the  general  form  of  the 
governing  two-dimensional  equations  was  discussed.  The  following  section 
represents  a  complete  list  of  the  necessary  equations  for  solving  the 
finite  difference  analysis  of  head  impact.  Each  grid  cone  utilizes 
similar  wave  propagation  relationships.  The  equations  given  below  are 
for  a  typical  elastic  material  grid  zone. 

The  nodal  location  identifiers  and  all  notations  are  described  in 
Figure  A-l  and  List  of  Symbols.  It  is  assumed  that  values  for  all  vari¬ 
ables  are  known  at  time  n.  The  system  of  equations  given  below  computes 
the  new  values  for  all  variables  at  time  n+1  for  the  1st  zone  called 
region  3  (also  defined  by  nodal  point  0).  After  completing  the  solution 
for  a  particular  zone,  the  next  zone  is  temporarily  identified  as  region 
3,  and  the  computations  are  repeated. 


Figure  Al.  Grid  Zone  Description  for  Finite 
Difference  Algorithm 
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A1 


The  conservation  of  linear  monentun  expressions  yield  the  follow¬ 
ing  acceleration  values  for  grid  zone  3  at  tine  n  (note:  the  stress 
values  are  the  sun  of  the  constitutive  stress  terns  plus  the  artificial 
viscous  stress  terns): 
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The  value  for  pa  is  obtained  from  the  value  of  the  regions  surround¬ 
ing  point  0  as 

PA  -  ^  [PiAi  ♦  P2A2  +P3A3  +  P4A4]  • 

The  velocity  of  point  0  is  advanced  in  time  by 
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Similarly,  the  new  location  of  point  0  is  given  by 
x"*1  -  X*  .  4(n*1/2  (UX)"*V2 
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The  stretching  deviators  contain  local  strain  rate  and  dilatation 
components .  The  final  difference  foim  for  these  deviators  at  time 

n+1/2  is  given  by  (Note:  ^d^'  is  eliminated) 

x-,  Xx  •,  ,  n+T  n  n+1  n. 
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Equations  A1,A2,A3  and  A6  require  density  and  area  values  that  are 
obtained  from  the  following  expressions: 
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The  assumed  energy  independent  nature  of  the  constitutive  equa¬ 
tions  yields  the  following  hydrodynamic  pressure  relationship 

.p5+1  «  Knn+*.  AB 

the  stress  deviator  amputations  require  the  rigid  body  rotation 
(W*2)  at  time  n+1/2  given  as 
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Now  using  Hooke's  Law  the  deviator  stress  components  can  be  computed  as 
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For  the  present  analysis  only  bulk  conponents  of  the  artificial 
viscosity  tensor  are  used.  Therefore,  the  necessary  terms  for  artifi¬ 
cial  stress  are  derived  as 
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Finally,  the  stress  conponents  at  time  n+1  are  conputed  by 
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The  remaining  stress  term  is  solved  as 
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